# Call the data cleaning script
source("./code/data_cleaning.R")Warning: package 'groundhog' was built under R version 4.4.1
Attached: 'Groundhog' (Version: 3.2.1)
Tips and troubleshooting: https://groundhogR.com
The following tracks the data analysis for our 2024 Undergraduate Student Research Award (USRA) project, “What accounts for the relationship between disgust and religiosity?” This is the master analysis document, so all analyses are documented here—with any other code called within a chunk below. The goals of this analysis are to (1) clean the data, (2) check for internal reliability of our measures, (3) generate descriptive statistics, (4) visualize the distribution of our variables, (5) test our hypotheses, and (6) conduct exploratory analyses as appropriate.
Individual differences in disgust sensitivity have repeatedly been shown to relate to religiosity, as predicted by the behavioral immune system model of disgust, leading to some scholars suggesting that religiosity is an evolved disease avoidance strategy. However, there are multiple mechanisms that aim to account for how religion contributes to disease avoidance, and these different mechanisms provide predictions for the psychological mediators of the relationship between disgust and religiosity. This study aims to test three such mediators, to better differentiate and understand why disgust sensitivity relates to religiosity. These mediators are represented by the following hypotheses.
Strong version of statistical hypothesis:
Weak version of statistical hypothesis:
Strong version of statistical hypothesis:
Weak version of statistical hypothesis:
Strong version of statistical hypothesis:
Weak version of statistical hypothesis:
We will test these hypotheses through mediation analysis, the specifics of which will be reviewed in the Testing Hypotheses section.
All data files are saved in .csv format in the data directory, and the steps that were taken to process these files are documented in data/data_preparation_notes.txt. To clean the data for analysis, we will call a pre-written script (code/data_cleaning.R) that returns a data frame called data. The script also saves this data frame as data/data_clean.csv.
# Call the data cleaning script
source("./code/data_cleaning.R")Warning: package 'groundhog' was built under R version 4.4.1
Attached: 'Groundhog' (Version: 3.2.1)
Tips and troubleshooting: https://groundhogR.com
Note: The Three Domains of Disgust Scale (TDDS) is generally administered with seven response options (Olatunji et al., 2012), only anchored with 0 (Not disgusting at all) and 6 (Extremely disgusting). Due to a survey coding error, our sample completed the TDDS with only six response options and the same anchors. Because of this, each TDDS item has a range from 0 to 5 (rather than 6), and raw scores for the full scale and subscales will be underestimated accordingly.
The data data frame contains the following variables (and column names).
Participant ID (ID):
Start Time (start_time) and Finish Time (finish_time):
Age (age):
Sex (sex):
Gender (gender):
The gender of the participant. Responses were standardized into six categories formulated based on the genders reported. To see which original responses were put into these bins, see the appropriate section of the code/data_cleaning.R script. The categories are as follows:
Ethnicity (ethnicity):
The ethnicity of the participant. Participants were provided with a list of potential ethnicities to choose from as well as the options of writing in their own or “Rather Not Say”. Those who reported multiple ethnicities were categorized as “Mixed”. There are 11 unique ethnicities for this sample:
Nationality (nationality):
Religious Affiliation (religious_affiliation):
Christian Affiliation (christian_affiliation):
Religious Importance (religious_importance):
The Centrality of Religiosity Scale (CRS):
The Three Domains of Disgust Scale (full, TDDS_f; pathogen, TDDS_p; sexual, TDDS_s, moral, TDDS_m):
The TDDS was developed by Tybur et al. (2009) in response to a lack of theoretical and empirical grounding for prior measures of disgust sensitivity. Based on evolutionary theory, they predicted that items which people find disgusting should form three factors, one for pathogen, sexual, and moral disgust. Across multiple studies, they demonstrate that a three-factor solution is most parsimonious and fits the data well, that the full scale has good concurrent validity, and that the subscales have good convergent and divergent validity. Subsequent investigations have confirmed the three factor structure of the scale, although they find poorer validity evidence for the moral disgust subscale (Olatunji et al., 2012). All three subscales will be used in exploratory analyses and bivariate correlations, but the pathogen disgust scale will be used to test our central three hypotheses and six specific statistical hypotheses.
Scores for the full scale and each subscale were calculated by summing item values.
The Revised Sociosexual Orientation Inventory (full, SOI_f; behavior, SOI_b; attitudes, SOI_a; desire, SOI_d):
The original Sociosexual Orientation Inventory (SOI) was a 7-item global measure of sociosexual orientation (Simpson & Gangestad, 1991), but this measure included behavioral, attitudinal, and (one) desire items together, with different response formats. This is problematic because each of these three components are differentiable both conceptually and empirically, and the different response formats often lead to an improper weighting of each of the three components contributing to the total score, resulting in poor reliability. Penke & Asendorpf (2008) developed the SOI-R to address these, and other, problems. The SOI-R is a 9-item measure with three items each tailored towards behaviors, attitudes, and desires, which they differentiated based on exploratory and confirmatory factor analysis as well as convergent and divergent validation. Each of the subscales show adequate reliability, along with the global score. As we are interested in operationalizing a monogamous mating strategy, it is participants intentions that we were most interested in measuring. In this case, one’s desire and behavior may, or may not, correspond to their intentions, which would best be operationalized as sociosexual attitudes. Therefore, we will use the attitudes subscale of the SOI-R to operationalize a monogamous mating strategy. This will be used to test our mediational hypothesis involving monogamous sexual strategies.
Scores for the full scale and each subscale were calculated by taking the arithmetic mean of item values.
The In-Group Preference Subscale of the Short Form of the Generalized Ethnocentrism Scale (SFGENE-7; full, GENE_f; preference, GENE_p; superiority, GENE_s):
The SFGENE-7 is simply a short form of the Generalized Ethnocentrism Scale (Neuliep, 2002; Neuliep & McCroskey, 1997). The GENE is composed of 22 items, 15 of which are scored, and out of all the measures that we have, the GENE would be by far the longest. By reducing the number of items to 7 with the SFGENE-7, we can allow for the recruitment of more participants. The SFGENE-7 has shown adequate reliability despite a loss of items, as well as evidence for convergent validity with constructs like tolerance and multicultural ideology. In addition, it has been broken down via factor analysis into two distinct subscales, the first three items comprising the in-group preference scale and the last four items comprosing the in-group superiority scale. Each of these subscales also have strong internal reliabilities, and the in-group preference scale shows stronger negative correlations with tolerance and multicultural identity than the in-group superiority scale. We will use use only the first three items to measure the in-group preference aspect of ethnocentrism in order to operationalize out-group avoidance, particularly because the items comprising the in-group superiority subscale are quite similar at face value to items in the conventionalism scale that we are using to measure the tendency towards adherence to tradition. The in-group preference subscale will be used to test the mediational hypothesis involving out-group avoidance.
Items were summed to provide the total, preference, and superiority scores.
Conventionalism Subscale of Aggression-Submission-Conventionalism Scale of Authoritarianism (ASC-C; CONV) and additional Traditionalism Items (TRAD; together, CONV_f):
TRAD items in the data frame. We will calculate Chronbach’s alpha and only accept these items in the measure if we keep an alpha above .8. If the alpha is below .8, we will remove items until it performs at this level.To quickly assess the internal reliability of our scales, we will calculate Chronbach’s alpha for each of our scales. We will start with the CONV_f scale, because we may remove items we added to it.
# Defining the packages to use to do the alpha analysis
pkg <- c("tidyverse", "psych")
# tidyverse if not attached from data cleaning script (for pipe and select function)
# psych for alpha (and later descriptives)
# Loading the groundhog package to install packages from a certain date
library(groundhog)
# Reading in the packages with the groundhog package for reproducibility
# Message suppressed in order to allow for rendering of quarto document
suppressMessages(groundhog.library(pkg = pkg, date = "2024-06-01"))
# Remove the pkg vector
rm(pkg)# Chronbach's Alpha for CONV_f (which includes all CONV and TRAD items)
alpha(x = data %>% select(starts_with("CONV."), starts_with("TRAD.")))
Reliability analysis
Call: alpha(x = data %>% select(starts_with("CONV."), starts_with("TRAD.")))
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.83 0.83 0.87 0.29 4.9 0.015 3 0.86 0.28
95% confidence boundaries
lower alpha upper
Feldt 0.8 0.83 0.86
Duhachek 0.8 0.83 0.86
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
CONV.1 0.82 0.82 0.86 0.30 4.6 0.015 0.034 0.29
CONV.2 0.80 0.80 0.84 0.27 4.0 0.017 0.031 0.25
CONV.3 0.82 0.82 0.85 0.29 4.6 0.015 0.029 0.30
CONV.4 0.81 0.81 0.85 0.28 4.4 0.016 0.034 0.28
CONV.5 0.82 0.81 0.85 0.28 4.4 0.016 0.034 0.28
CONV.6 0.82 0.82 0.85 0.29 4.5 0.016 0.028 0.29
TRAD.1 0.84 0.84 0.87 0.33 5.4 0.014 0.023 0.31
TRAD.2 0.82 0.82 0.85 0.29 4.5 0.016 0.030 0.28
TRAD.3 0.82 0.82 0.86 0.29 4.5 0.015 0.034 0.28
TRAD.4 0.82 0.81 0.85 0.28 4.4 0.016 0.031 0.29
TRAD.5 0.81 0.80 0.84 0.27 4.1 0.017 0.028 0.26
TRAD.6 0.81 0.81 0.85 0.28 4.3 0.016 0.030 0.28
Item statistics
n raw.r std.r r.cor r.drop mean sd
CONV.1 289 0.53 0.53 0.47 0.42 2.3 1.4
CONV.2 289 0.76 0.75 0.74 0.68 3.4 1.5
CONV.3 289 0.54 0.55 0.52 0.43 3.4 1.5
CONV.4 289 0.64 0.63 0.58 0.53 2.7 1.5
CONV.5 289 0.63 0.63 0.59 0.54 2.2 1.5
CONV.6 289 0.57 0.58 0.55 0.47 3.9 1.3
TRAD.1 289 0.27 0.27 0.17 0.14 1.5 1.4
TRAD.2 289 0.57 0.57 0.53 0.47 3.9 1.4
TRAD.3 289 0.57 0.57 0.51 0.46 2.0 1.5
TRAD.4 289 0.63 0.63 0.60 0.53 3.0 1.5
TRAD.5 289 0.73 0.73 0.73 0.66 3.7 1.4
TRAD.6 289 0.63 0.64 0.61 0.54 3.7 1.4
Non missing response frequency for each item
0 1 2 3 4 5 6 miss
CONV.1 0.13 0.17 0.27 0.25 0.12 0.03 0.02 0
CONV.2 0.04 0.07 0.14 0.28 0.24 0.13 0.10 0
CONV.3 0.04 0.06 0.12 0.29 0.25 0.15 0.09 0
CONV.4 0.07 0.13 0.27 0.26 0.13 0.08 0.06 0
CONV.5 0.13 0.18 0.31 0.21 0.10 0.04 0.03 0
CONV.6 0.01 0.03 0.09 0.25 0.30 0.19 0.12 0
TRAD.1 0.31 0.25 0.22 0.15 0.03 0.02 0.01 0
TRAD.2 0.02 0.04 0.10 0.20 0.31 0.17 0.15 0
TRAD.3 0.19 0.23 0.25 0.16 0.11 0.04 0.01 0
TRAD.4 0.06 0.09 0.23 0.25 0.23 0.09 0.06 0
TRAD.5 0.03 0.05 0.09 0.26 0.32 0.13 0.12 0
TRAD.6 0.02 0.07 0.08 0.25 0.30 0.20 0.08 0
The raw alpha value (.831; standardized .830) is acceptable here, so we will use our full CONV_f variable to operationalize traditionalism.
To compare these values to those for the original conventionalism measure, I will calculate its’ alpha as well.
# Chronbach's Alpha for CONV (only CONV items)
alpha(x = data %>% select(starts_with("CONV.")))
Reliability analysis
Call: alpha(x = data %>% select(starts_with("CONV.")))
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.74 0.74 0.76 0.33 2.9 0.024 3 0.97 0.38
95% confidence boundaries
lower alpha upper
Feldt 0.7 0.74 0.79
Duhachek 0.7 0.74 0.79
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
CONV.1 0.73 0.73 0.75 0.36 2.8 0.025 0.028 0.40
CONV.2 0.67 0.67 0.70 0.29 2.1 0.031 0.038 0.21
CONV.3 0.72 0.72 0.70 0.34 2.6 0.025 0.018 0.39
CONV.4 0.70 0.70 0.71 0.32 2.4 0.028 0.031 0.36
CONV.5 0.69 0.69 0.70 0.31 2.3 0.029 0.030 0.35
CONV.6 0.72 0.71 0.70 0.33 2.5 0.026 0.020 0.39
Item statistics
n raw.r std.r r.cor r.drop mean sd
CONV.1 289 0.59 0.59 0.45 0.39 2.3 1.4
CONV.2 289 0.75 0.75 0.67 0.59 3.4 1.5
CONV.3 289 0.62 0.63 0.56 0.43 3.4 1.5
CONV.4 289 0.68 0.67 0.58 0.50 2.7 1.5
CONV.5 289 0.70 0.70 0.62 0.53 2.2 1.5
CONV.6 289 0.63 0.64 0.58 0.45 3.9 1.3
Non missing response frequency for each item
0 1 2 3 4 5 6 miss
CONV.1 0.13 0.17 0.27 0.25 0.12 0.03 0.02 0
CONV.2 0.04 0.07 0.14 0.28 0.24 0.13 0.10 0
CONV.3 0.04 0.06 0.12 0.29 0.25 0.15 0.09 0
CONV.4 0.07 0.13 0.27 0.26 0.13 0.08 0.06 0
CONV.5 0.13 0.18 0.31 0.21 0.10 0.04 0.03 0
CONV.6 0.01 0.03 0.09 0.25 0.30 0.19 0.12 0
CONV_f variable has stronger internal reliability than the CONV items alone.# Chronbach's Alpha for Centrality of Religiosity Items
alpha(x = data %>% select(starts_with("CRS.")))
Reliability analysis
Call: alpha(x = data %>% select(starts_with("CRS.")))
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.92 0.92 0.91 0.69 11 0.0069 2.6 1.3 0.72
95% confidence boundaries
lower alpha upper
Feldt 0.90 0.92 0.93
Duhachek 0.91 0.92 0.93
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
CRS.1 0.93 0.93 0.92 0.78 13.8 0.0064 0.002 0.80
CRS.2 0.89 0.89 0.88 0.68 8.4 0.0094 0.012 0.66
CRS.3 0.90 0.90 0.88 0.69 8.8 0.0090 0.016 0.70
CRS.4 0.89 0.89 0.86 0.66 7.7 0.0104 0.011 0.65
CRS.5 0.89 0.89 0.87 0.67 8.1 0.0096 0.015 0.65
Item statistics
n raw.r std.r r.cor r.drop mean sd
CRS.1 289 0.73 0.76 0.65 0.63 2.6 1.1
CRS.2 289 0.90 0.89 0.87 0.83 3.0 1.5
CRS.3 289 0.88 0.88 0.84 0.81 2.3 1.4
CRS.4 289 0.93 0.92 0.91 0.87 2.5 1.7
CRS.5 289 0.90 0.90 0.88 0.85 2.5 1.4
Non missing response frequency for each item
1 2 3 4 5 miss
CRS.1 0.14 0.37 0.28 0.13 0.08 0
CRS.2 0.18 0.29 0.15 0.11 0.27 0
CRS.3 0.42 0.22 0.13 0.09 0.14 0
CRS.4 0.46 0.15 0.08 0.03 0.28 0
CRS.5 0.33 0.26 0.16 0.11 0.14 0
# Chronbach's Alpha for Pathogen Disgust Items
alpha(x = data %>% select(., c("TDDS.12", "TDDS.15", "TDDS.9", "TDDS.3", "TDDS.21", "TDDS.18", "TDDS.6")))
Reliability analysis
Call: alpha(x = data %>% select(., c("TDDS.12", "TDDS.15", "TDDS.9",
"TDDS.3", "TDDS.21", "TDDS.18", "TDDS.6")))
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.77 0.77 0.76 0.33 3.4 0.021 3 0.91 0.34
95% confidence boundaries
lower alpha upper
Feldt 0.73 0.77 0.81
Duhachek 0.73 0.77 0.81
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
TDDS.12 0.74 0.74 0.72 0.33 2.9 0.023 0.0094 0.34
TDDS.15 0.72 0.72 0.70 0.30 2.6 0.026 0.0063 0.28
TDDS.9 0.72 0.73 0.70 0.31 2.7 0.025 0.0064 0.31
TDDS.3 0.75 0.75 0.74 0.34 3.1 0.023 0.0091 0.34
TDDS.21 0.75 0.75 0.74 0.33 3.0 0.023 0.0114 0.35
TDDS.18 0.75 0.75 0.74 0.34 3.1 0.022 0.0095 0.35
TDDS.6 0.75 0.76 0.74 0.34 3.1 0.023 0.0076 0.34
Item statistics
n raw.r std.r r.cor r.drop mean sd
TDDS.12 289 0.65 0.65 0.57 0.49 3.4 1.4
TDDS.15 289 0.73 0.74 0.70 0.61 3.3 1.3
TDDS.9 289 0.71 0.71 0.67 0.57 2.5 1.4
TDDS.3 289 0.58 0.61 0.51 0.44 4.0 1.1
TDDS.21 289 0.65 0.63 0.52 0.47 3.0 1.5
TDDS.18 289 0.63 0.61 0.51 0.44 3.1 1.5
TDDS.6 289 0.61 0.60 0.50 0.44 1.9 1.4
Non missing response frequency for each item
0 1 2 3 4 5 miss
TDDS.12 0.02 0.09 0.15 0.21 0.25 0.28 0
TDDS.15 0.01 0.11 0.15 0.26 0.24 0.24 0
TDDS.9 0.07 0.20 0.28 0.21 0.15 0.10 0
TDDS.3 0.00 0.03 0.09 0.17 0.23 0.48 0
TDDS.21 0.06 0.15 0.18 0.22 0.15 0.24 0
TDDS.18 0.04 0.14 0.18 0.19 0.17 0.27 0
TDDS.6 0.17 0.24 0.26 0.18 0.08 0.06 0
TDDS_p) shows adequate internal reliability here (alpha = .769; standardized = .772).# Chronbach's Alpha for the Attitudes Subscale of the SOI-R
alpha(x = data %>% select(., c("SOI.4", "SOI.5", "SOI.6")))
Reliability analysis
Call: alpha(x = data %>% select(., c("SOI.4", "SOI.5", "SOI.6")))
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.81 0.81 0.74 0.58 4.2 0.02 5.2 2.4 0.6
95% confidence boundaries
lower alpha upper
Feldt 0.76 0.81 0.84
Duhachek 0.77 0.81 0.84
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
SOI.4 0.76 0.76 0.62 0.62 3.2 0.028 NA 0.62
SOI.5 0.69 0.69 0.52 0.52 2.2 0.037 NA 0.52
SOI.6 0.75 0.75 0.60 0.60 3.0 0.029 NA 0.60
Item statistics
n raw.r std.r r.cor r.drop mean sd
SOI.4 289 0.83 0.83 0.70 0.63 5.9 2.7
SOI.5 289 0.88 0.87 0.78 0.70 4.3 2.9
SOI.6 289 0.84 0.84 0.71 0.64 5.3 2.8
Non missing response frequency for each item
1 2 3 4 5 6 7 8 9 miss
SOI.4 0.11 0.04 0.08 0.06 0.11 0.09 0.15 0.11 0.25 0
SOI.5 0.29 0.11 0.10 0.04 0.10 0.09 0.09 0.06 0.13 0
SOI.6 0.17 0.04 0.11 0.05 0.14 0.09 0.08 0.12 0.19 0
SOI_a) seems to have good internal consistency here (alpha = .806; standardized = .806).# Chronbach's Alpha for the Preferences Subscale of the GENE
alpha(x = data %>% select(., c("GENE.1", "GENE.2", "GENE.3")))
Reliability analysis
Call: alpha(x = data %>% select(., c("GENE.1", "GENE.2", "GENE.3")))
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.52 0.53 0.55 0.27 1.1 0.05 2.3 0.68 0.095
95% confidence boundaries
lower alpha upper
Feldt 0.42 0.52 0.61
Duhachek 0.43 0.52 0.62
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
GENE.1 0.801 0.804 0.672 0.672 4.09 0.023 NA 0.672
GENE.2 0.093 0.093 0.049 0.049 0.10 0.107 NA 0.049
GENE.3 0.174 0.174 0.095 0.095 0.21 0.097 NA 0.095
Item statistics
n raw.r std.r r.cor r.drop mean sd
GENE.1 289 0.54 0.53 0.096 0.078 3.2 0.95
GENE.2 289 0.81 0.82 0.766 0.536 1.7 0.89
GENE.3 289 0.80 0.80 0.734 0.472 1.9 0.99
Non missing response frequency for each item
1 2 3 4 5 miss
GENE.1 0.05 0.12 0.47 0.27 0.09 0
GENE.2 0.51 0.32 0.12 0.03 0.01 0
GENE.3 0.47 0.29 0.18 0.06 0.01 0
GENE_p_2.# Creating the GENE_p_2 variable as the sum of items 2 and 3
data <- data %>%
mutate(
GENE_p_2 = rowSums(select(., c(GENE.2, GENE.3)))
)First, we will see how many participants we have and how long it took participants to complete the survey.
# Calculate the number of participants in the sample
nrow(data)[1] 289
# Calculate descriptive statistics for the difference between the finish and start time (in minutes)
describe(as.numeric(difftime(data$finish_time, data$start_time, units = "mins"))) vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 289 14.39 58.1 8 8.94 2.97 3 821 818 12.25 154.08 3.42
In total, we have N = 289 participants with complete data.
Although the mean time taken is around 14.39, the SD is huge, and the maximum is a massive 821 minutes (over 13 hours). The trimmed mean of 8.94 is therefore a better estimate.
Next, we will describe the sample in terms of the categorical variables. The following generates frequency tables for all categorical variables.
# Create a frequency table for each factor variable
lapply(data[sapply(data, is.factor)], table)$sex
Male Female
146 143
$gender
Agender Female Gender Fluid Gender Queer Male Non-binary
2 135 1 1 144 4
$ethnicity
African Black or African American
52 4
Caribbean East Asian
0 14
Latino or Hispanic Middle Eastern
62 4
Mixed Native American or Alaskan Native
11 0
South Asian White or Causasian
8 122
White or Sapharic Jew Black British
6 1
White Mexican Romani or Traveller
0 0
South East Asian Rather Not Say
2 3
$nationality
Algeria Australia Austria Belgium Brazil
1 7 1 1 3
Canada Chile China Columbia Czech Republic
30 13 1 2 5
Egypt Eritrea Finland France Germany
1 1 2 3 4
Ghana Greece Hungary India Ireland
1 5 3 1 1
Israel Italy Japan Kenya Mexico
2 5 1 6 39
New Zealand Poland Portugal Romania Saudi Arabia
6 14 25 1 1
Slovakia South Africa Spain Sweden Tunisia
1 46 7 5 1
Turkey United Kingdom United States Venezuela Vietnam
2 21 17 1 1
Zimbabwe
1
$religious_affiliation
Agnostic Anti-religious Atheist Buddhist Christian
6 1 2 2 123
Deist Hindu Muslim None Pagan
1 1 11 139 1
Spiritual Spiritualist
1 1
$christian_affiliation
Anglican Apostolic Baptist
1 1 11
Catholic Congregationalist Evangelical
68 1 4
Jehovah’s Witness Lutheran Methodist
1 6 9
Non-religious Christian None Orthodox
1 4 4
Pentecostal Presbyterian Protestant
1 2 8
Seventh Day Adventist
1
The output is likely difficult to read, so I will break it down here.
Sex:
Gender:
Male: 144 (49.8%)
Female: 135 (46.7%)
Non-Binary: 4 (1.4%)
Agender: 2 (.7%)
Gender Fluid: 1 (.3%)
Gender Queer: 1 (.3%)
Ethnicity:
Nationality:
Religious Affiliation:
Christian Affiliation:
Now we will describe our numeric variables starting with our only demographic numeric variable, age.
# Generate descriptive statistics for age
describe(data$age) vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 289 30.85 10.85 28 29.18 7.41 18 86 68 1.74 3.84 0.64
Now we will describe the important scales for our analysis.
# Describe important scales in our analysis
describe(data %>% select("religious_importance", "CRS", "TDDS_f", "TDDS_p", "SOI_f", "SOI_a", "GENE_p", "GENE_p_2", "CONV_f")) vars n mean sd median trimmed mad min max range
religious_importance 1 289 30.70 35.63 12.00 26.08 17.79 0 100 100
CRS 2 289 2.59 1.26 2.20 2.50 1.19 1 5 4
TDDS_f 3 289 60.12 16.05 59.00 59.82 16.31 12 100 88
TDDS_p 4 289 21.26 6.35 22.00 21.30 5.93 4 35 31
SOI_f 5 289 3.58 1.62 3.22 3.48 1.81 1 8 7
SOI_a 6 289 5.16 2.40 5.00 5.19 2.47 1 9 8
GENE_p 7 289 6.81 2.03 7.00 6.72 1.48 3 15 12
GENE_p_2 8 289 3.57 1.71 3.00 3.34 1.48 2 10 8
CONV_f 9 289 35.55 10.30 36.00 35.78 10.38 0 65 65
skew kurtosis se
religious_importance 0.87 -0.83 2.10
CRS 0.57 -1.06 0.07
TDDS_f 0.17 -0.29 0.94
TDDS_p -0.07 -0.49 0.37
SOI_f 0.53 -0.59 0.10
SOI_a -0.03 -1.07 0.14
GENE_p 0.50 0.36 0.12
GENE_p_2 1.02 0.64 0.10
CONV_f -0.24 0.44 0.61
religious_importance and the GENE_p_2, there does not seem to be a problem here. Variables will be inspected more closely when testing assumptions for statistical tests.Before moving on, we will calculate correlations between each of these variables.
# Calculate correlations for numerical variables
corr.test(data %>% select("age", "religious_importance", "CRS", "TDDS_f", "TDDS_p", "TDDS_s", "SOI_f", "SOI_a", "GENE_p", "GENE_p_2", "CONV_f"))Call:corr.test(x = data %>% select("age", "religious_importance",
"CRS", "TDDS_f", "TDDS_p", "TDDS_s", "SOI_f", "SOI_a", "GENE_p",
"GENE_p_2", "CONV_f"))
Correlation matrix
age religious_importance CRS TDDS_f TDDS_p TDDS_s
age 1.00 -0.07 -0.09 0.02 -0.08 -0.09
religious_importance -0.07 1.00 0.92 0.46 0.22 0.52
CRS -0.09 0.92 1.00 0.47 0.24 0.53
TDDS_f 0.02 0.46 0.47 1.00 0.76 0.82
TDDS_p -0.08 0.22 0.24 0.76 1.00 0.47
TDDS_s -0.09 0.52 0.53 0.82 0.47 1.00
SOI_f 0.12 -0.28 -0.27 -0.38 -0.16 -0.54
SOI_a 0.14 -0.39 -0.41 -0.45 -0.22 -0.59
GENE_p -0.06 0.07 0.07 0.02 0.01 0.07
GENE_p_2 0.03 0.02 0.00 -0.02 -0.03 0.04
CONV_f 0.08 0.33 0.31 0.12 0.01 0.13
SOI_f SOI_a GENE_p GENE_p_2 CONV_f
age 0.12 0.14 -0.06 0.03 0.08
religious_importance -0.28 -0.39 0.07 0.02 0.33
CRS -0.27 -0.41 0.07 0.00 0.31
TDDS_f -0.38 -0.45 0.02 -0.02 0.12
TDDS_p -0.16 -0.22 0.01 -0.03 0.01
TDDS_s -0.54 -0.59 0.07 0.04 0.13
SOI_f 1.00 0.88 -0.11 -0.10 -0.10
SOI_a 0.88 1.00 -0.15 -0.11 -0.15
GENE_p -0.11 -0.15 1.00 0.88 0.08
GENE_p_2 -0.10 -0.11 0.88 1.00 0.03
CONV_f -0.10 -0.15 0.08 0.03 1.00
Sample Size
[1] 289
Probability values (Entries above the diagonal are adjusted for multiple tests.)
age religious_importance CRS TDDS_f TDDS_p TDDS_s SOI_f
age 0.00 1.00 1.00 1.00 1.00 1.00 1.00
religious_importance 0.21 0.00 0.00 0.00 0.01 0.00 0.00
CRS 0.15 0.00 0.00 0.00 0.00 0.00 0.00
TDDS_f 0.79 0.00 0.00 0.00 0.00 0.00 0.00
TDDS_p 0.15 0.00 0.00 0.00 0.00 0.00 0.21
TDDS_s 0.14 0.00 0.00 0.00 0.00 0.00 0.00
SOI_f 0.05 0.00 0.00 0.00 0.01 0.00 0.00
SOI_a 0.01 0.00 0.00 0.00 0.00 0.00 0.00
GENE_p 0.35 0.20 0.25 0.75 0.81 0.21 0.06
GENE_p_2 0.66 0.79 1.00 0.72 0.67 0.55 0.09
CONV_f 0.17 0.00 0.00 0.04 0.93 0.03 0.10
SOI_a GENE_p GENE_p_2 CONV_f
age 0.42 1.00 1.00 1.00
religious_importance 0.00 1.00 1.00 0.00
CRS 0.00 1.00 1.00 0.00
TDDS_f 0.00 1.00 1.00 1.00
TDDS_p 0.01 1.00 1.00 1.00
TDDS_s 0.00 1.00 1.00 0.74
SOI_f 0.00 1.00 1.00 1.00
SOI_a 0.00 0.29 1.00 0.32
GENE_p 0.01 0.00 0.00 1.00
GENE_p_2 0.06 0.00 0.00 1.00
CONV_f 0.01 0.16 0.62 0.00
To see confidence intervals of the correlations, print with the short=FALSE option
For each correlation the sample size is N = 289.
Noteworthy for this analysis, pathogen disgust is positively and significantly related to both religious_importance (r = .22, p < .01) and the CRS (r = .24, p < .01), but sexual disgust is more strongly related to both religious_importance (r = .52, p < .01) and the CRS (r = .53, p < .01). This replicates previous work and provides initial support for the sexual strategies hypothesis.
In addition, there is a very strong correlation between religious_importance and the CRS (r = .92, p < .01).
However, neither the GENE_p nor the GENE_p_2 are significantly related to either religious_importance (r = .07, p = .20; r = .02, p = .02; respectively) or the CRS (r = .07, p = .25; r = .00, p = .00; respectively), suggesting that it could not mediate the relationship between disgust sensitivity and religiosity.
Unlike the GENE variables, CONV_f is significantly related to both religious_importance (r = .33, p < .01) and the CRS (r = .31, p < .01).
Next, the only potential mediator that is significantly related to pathogen disgust is the SOI_a:
SOI_a: r = -.22, p = .01
GENE_p: r = .01, p = 1
GENE_p_2: r = -.03, p = 1CONV_f: r = .01, p = 1
This suggests that the only plausible mediator (here, that is) of the relationship between disgust and religiosity could be monogamous mating strategies.
Now we will visualization the distribution of our variables by creating bar plots and histograms (for categorical and numerical variables, respectively).
# Loop through each categorical variable in the data frame and create a bar chart
for (var in names(data[, sapply(data, is.factor)])) {
p <- ggplot(data[, sapply(data, is.factor)], aes(x = !!sym(var))) +
geom_bar() +
theme_minimal() +
labs(title = paste("Bar Chart of", var), x = var, y = "Count") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8)) # Adjust x axis labels to make them readable (not all over eachother)
# Print each bar chart
print(p)
# Remove the plot object so it doesn't clutter the environment
rm(p)
}# Remove the var object
rm(var)# Loop through numeric variable identified to produce a histogram
for (var in names(data %>% select("age", "religious_importance", "CRS", "TDDS_f", "TDDS_p", "SOI_f", "SOI_a", "GENE_p", "GENE_p_2", "CONV_f"))) {
p <- ggplot(data %>% select("age", "religious_importance", "CRS", "TDDS_f", "TDDS_p", "SOI_f", "SOI_a", "GENE_p", "GENE_p_2", "CONV_f"), aes(x = .data[[var]])) +
geom_histogram(bins = 10, fill = "steelblue", color = "black") +
theme_minimal() +
labs(title = paste("Histogram of", var), x = var, y = "Frequency")
# Print each histogram
print(p)
# Remove the p object so it doesn't clutter up the environment
rm(p)
}# Remove the var object from the environment
rm(var)As above, the following are our hypotheses.
Strong version of statistical hypothesis:
Weak version of statistical hypothesis:
Strong version of statistical hypothesis:
Weak version of statistical hypothesis:
Strong version of statistical hypothesis:
Weak version of statistical hypothesis:
The three hypotheses are each broken down into two statistical hypotheses, which represent their strong and weak version. In the strong version there is full mediation, and in the weak version there is only partial mediation. We have defined and preregistered inference criteria for both of these versions of each hypothesis.
Strong Version Inference Criterion:
For each model (each hypothesis) we will use the “K” method outlined by Beribisky et al. (2020) to infer full mediation. That is, if the proportion of variance explained by the indirect effect compared to the total effect is 80% or higher, we will infer full mediation. This will provide support for the strong version of the respective hypothesis.
Weak Version Inference Criterion:
Where there is not full mediation, for each model (each hypothesis) if the 95% confidence interval for the indirect effect of parasite disgust sensitivity on religiosity through the mediator does not contain zero, then we will infer partial mediation. This will provide support for the weak version of the respective hypothesis.
Before testing our hypotheses, we will standardize each variable. We will also reverse the scores of SOI_a before standardizing it, because higher scores currently represent more unrestricted sociosexual attitudes, whereas we would like higher scores to represent more restricted sociosexual attitudes (as a proxy for a monogamous mating strategy).
# Reversing the scoring for the SOI_a (saving it as SOI_a_r) and standardizing all variables
data <- data %>%
mutate(
SOI_a_r = (max(data$SOI_a) + 1) - SOI_a,
SOI_a_r_z = scale(SOI_a_r),
CRS_z = scale(CRS),
religious_importance_z = scale(religious_importance),
GENE_p_z = scale(GENE_p),
GENE_p_2_z = scale(GENE_p_2),
CONV_f_z = scale(CONV_f),
TDDS_p_z = scale(TDDS_p)
)
# Ensure that the variables are plain numeric vectors
data$TDDS_p_z <- as.numeric(data$TDDS_p_z)
data$CONV_f_z <- as.numeric(data$CONV_f_z)
data$CRS_z <- as.numeric(data$CRS_z)
data$SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
data$religious_importance_z <- as.numeric(data$religious_importance_z)
data$GENE_p_z <- as.numeric(data$GENE_p_z)
data$GENE_p_2_z <- as.numeric(data$GENE_p_2_z)
# Because the boot function was throwing a fit about this earlier
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)We will also load the packages involved in the mediation analysis.
# Defining the packages to use to do the mediation analysis
pkg <- c("car", "mediation")
# car for testing for multicollinearity and bootstrapping
# mediation for testing the indirect effect with bootstrapped CIs and accessing
# Reading in the packages with the groundhog package for reproducibility
# Message suppressed in order to allow for rendering of quarto document
suppressMessages(groundhog.library(pkg = pkg, date = "2024-06-01"))
# Dropping the pkg vector
rm(pkg)Before moving on, we will first test the baseline model where we regress the dependent variable (CRS) on the independent variable (TDDS_p), because each mediation analysis requires that there is a significant positive relationship between these variables. We will also test each of the assumptions.
# Fitting the model
baseline_model <- lm(CRS_z ~ TDDS_p_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(baseline_model, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(baseline_model))# QQ-plot for normality of residuals
qqnorm(residuals(baseline_model))
qqline(residuals(baseline_model), col = "red")The histogram of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just around -1 SD and around 2 SD. The two clusters of similar errors in prediction may be due to sex differences in disgust and religiosity. In exploratory analyses, we will add sex as a control variable.
Because our inferential criteria do not require parametric tests, we should not be in too much trouble with a non-normality of residuals. We were relying on a parametric t-test to assess whether pathogen disgust is related to religiosity in the first place, however. Because we have failed the assumption of normality of residuals, we will use a bootstrapped test to make our inference here as well.
Multicollinearity: Not applicable becuase only one predictor
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(baseline_model, col = "black", lwd = 1)# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_baseline_model <- Boot(baseline_model, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_baseline_model) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -8.5501e-16 0.0015289 0.057618 0.0036089
TDDS_p_z 2.4299e-01 0.0010657 0.059478 0.2411963
confint(boot_baseline_model) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1197906 0.1067505
TDDS_p_z 0.1345033 0.3605490
# Removing the CRS_z and TDDS_p_z vectors in the environment that the Boot package required outside of the data.frame
rm(CRS_z, TDDS_p_z)We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and religiosity.
Based on the estimates, the standardized coefficient for pathogen disgust is around .24. This aligns with estimates from a meta-analysis of this relationship (Yu et al., 2022).
This shows that the IV effects the DV, as in the second regression equation for testing mediation in Baron & Kenny (1986). To establish mediation, we will still need to establish that the IV affects the mediator and the mediator affects the DV. Before we test the indirect effect and proportion of variance explained by each effect for each mediation model, we will run these two models.
As mentioned above, to test Hypothesis 1, we must first establish that pathogen disgust affects adherence to traditional practices (i.e., CONV_f_z) and that adherence to traditional practices influences religiosity.
# Fitting the models
H1_M1 <- lm(CONV_f_z ~ TDDS_p_z, data = data)
H1_M2 <- lm(CRS_z ~ CONV_f_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_M1, which = 1) # For model 1plot(H1_M2, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_M1)) # For model 1hist(residuals(H1_M2)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H1_M1))
qqline(residuals(H1_M1), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H1_M2))
qqline(residuals(H1_M2), col = "red")Model 1:
Model 2:
The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$CONV_f_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Traditionalism",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H1_M1, col = "black", lwd = 1)# Summarizing the model
summary(H1_M1)
Call:
lm(formula = CONV_f_z ~ TDDS_p_z, data = data)
Residuals:
Min 1Q Median 3Q Max
-3.4517 -0.6429 0.0402 0.6287 2.8524
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.379e-17 5.893e-02 0.000 1.000
TDDS_p_z 5.260e-03 5.903e-02 0.089 0.929
Residual standard error: 1.002 on 287 degrees of freedom
Multiple R-squared: 2.766e-05, Adjusted R-squared: -0.003457
F-statistic: 0.00794 on 1 and 287 DF, p-value: 0.9291
We can see from the Scatter Plot that there is no relationship between pathogen disgust and traditionalism to speak of (B = .005, t(287) = .089, p = .929), just as with the correlation above.
# Create a scatterplot with the regression line to visualize
plot(data$CONV_f_z, data$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Traditionalism",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H1_M2, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
CRS_z <- as.numeric(data$CRS_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M2 <- Boot(H1_M2, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_M2) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -8.5089e-16 0.0026934 0.05539 0.0029616
CONV_f_z 3.1443e-01 0.0036837 0.05224 0.3164595
confint(boot_H1_M2) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1128222 0.1044219
CONV_f_z 0.2036523 0.4178083
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z)We can see from the Scatter Plot that there is a positive relationship between traditionalism and religiosity.
Based on the estimates, the standardized coefficient for traditionalism is around .31.
The fact that pathogen disgust does not predict traditionalism indicates that traditionalism could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H1_full <- lm(CRS_z ~ TDDS_p_z + CONV_f_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_full, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_full))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H1_full))
qqline(residuals(H1_full), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H1_full)TDDS_p_z CONV_f_z
1.000028 1.000028
# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_full <- Boot(H1_full, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_full) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -8.3817e-16 0.00170100 0.053887 0.001676
TDDS_p_z 2.4134e-01 -0.00064514 0.056527 0.240610
CONV_f_z 3.1316e-01 0.00297386 0.048923 0.316369
confint(boot_H1_full) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1106346 0.1026964
TDDS_p_z 0.1328538 0.3560156
CONV_f_z 0.2100775 0.4057880
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H1_mediate_results <- mediate(model.m = H1_M1,
model.y = H1_full,
treat = "TDDS_p_z",
mediator = "CONV_f_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results for the mediation analysis
summary(H1_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.00165 -0.03393 0.04 0.96
ADE 0.24134 0.12018 0.35 <2e-16 ***
Total Effect 0.24299 0.13035 0.35 <2e-16 ***
Prop. Mediated 0.00678 -0.17652 0.19 0.96
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs not containing zero, both pathogen disgust and traditionalism are significantly positively related to religiosity.
The coefficients for pathogen disgust (.24) and traditionalism (.31) are virtually identical in this model to the models above where they are stand-alone predictors. This suggests that the unique effects of pathogen disgust and traditionalism on religiosity over-and-above each other are completely independent effects.
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through traditionalism (B = .001, BCI = [-.035, .04]).
This result is inconsistent with both the weak and strong version version of Hypothesis 1. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by traditionalism.
As mentioned above, to test Hypothesis 2, we must first establish that pathogen disgust affects out-group avoidance (i.e., GENE_p_z) and that out-group avoidance influences religiosity.
# Fitting the models
H2_M1 <- lm(GENE_p_z ~ TDDS_p_z, data = data)
H2_M2 <- lm(CRS_z ~ GENE_p_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_M1, which = 1) # For model 1plot(H2_M2, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_M1)) # For model 1hist(residuals(H2_M2)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H2_M1))
qqline(residuals(H2_M1), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H2_M2))
qqline(residuals(H2_M2), col = "red")Model 1:
Model 2:
The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 or 1.5 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$GENE_p_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Out-Group Avoidance (GENE_p)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M1, col = "black", lwd = 1)# Summarizing the model
summary(H2_M1)
Call:
lm(formula = GENE_p_z ~ TDDS_p_z, data = data)
Residuals:
Min 1Q Median 3Q Max
-1.9107 -0.8800 0.0804 0.5946 4.0239
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.976e-16 5.892e-02 0.000 1.000
TDDS_p_z 1.440e-02 5.902e-02 0.244 0.807
Residual standard error: 1.002 on 287 degrees of freedom
Multiple R-squared: 0.0002075, Adjusted R-squared: -0.003276
F-statistic: 0.05956 on 1 and 287 DF, p-value: 0.8074
We can see from the Scatter Plot that there is no relationship between pathogen disgust and out-group avoidance to speak of (B = .014, t(287) = .244, p = .807), just as with the correlation above.
# Create a scatterplot with the regression line to visualize
plot(data$GENE_p_z, data$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Out-Group Avoidance (GENE_p)",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M2, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
CRS_z <- as.numeric(data$CRS_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2 <- Boot(H2_M2, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_M2) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -8.8123e-16 0.00268445 0.059258 0.0010875
GENE_p_z 6.7778e-02 0.00038975 0.059761 0.0678462
confint(boot_H2_M2) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.11359766 0.1204258
GENE_p_z -0.04811381 0.1832814
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z)We can see from the Scatter Plot that there is a positive trend in the relationship between out-group avoidance and religiosity.
Based on the estimates, the standardized coefficient for out-group avoidance is around .068.
The fact that pathogen disgust does not predict out-group avoidance and, in turn, out-group avoidance does not predict religiosity indicates that out-group avoidance could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H2_full <- lm(CRS_z ~ TDDS_p_z + GENE_p_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_full, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_full))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H2_full))
qqline(residuals(H2_full), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H2_full)TDDS_p_z GENE_p_z
1.000208 1.000208
# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_full <- Boot(H2_full, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_full) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -8.6771e-16 0.00171383 0.057648 0.0032765
TDDS_p_z 2.4206e-01 0.00079979 0.058930 0.2414062
GENE_p_z 6.4291e-02 0.00032819 0.056867 0.0657697
confint(boot_H2_full) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.11852205 0.1083173
TDDS_p_z 0.13062951 0.3574411
GENE_p_z -0.05479052 0.1697520
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_mediate_results <- mediate(model.m = H2_M1,
model.y = H2_full,
treat = "TDDS_p_z",
mediator = "GENE_p_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results of the mediation analysis
summary(H2_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.000926 -0.008506 0.01 0.85
ADE 0.242061 0.130624 0.35 <2e-16 ***
Total Effect 0.242987 0.130352 0.35 <2e-16 ***
Prop. Mediated 0.003811 -0.046003 0.06 0.85
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs, pathogen disgust, but not out-group avoidance, is significantly positively related to religiosity.
The coefficients for pathogen disgust (.24) and out-group avoidance (.06) are virtually identical in this model to the models above where they are stand-alone predictors. This suggests that the unique effect of pathogen disgust (and not out-group avoidance) on religiosity over-and-above each other are completely independent.
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through out-group avoidance (B = .001, BCI = [-.009, .01]).
This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by out-group avoidance.
As mentioned above, to test Hypothesis 3, we must first establish that pathogen disgust affects one’s mating strategy (i.e., make it more monogamous; as measured by more restricted sociosexuality; SOI_a_r_z) and that a monogamous mating strategy influences religiosity.
# Fitting the models
H3_M1 <- lm(SOI_a_r_z ~ TDDS_p_z, data = data)
H3_M2 <- lm(CRS_z ~ SOI_a_r_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_M1, which = 1) # For model 1plot(H3_M2, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_M1)) # For model 1hist(residuals(H3_M2)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H3_M1))
qqline(residuals(H3_M1), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H3_M2))
qqline(residuals(H3_M2), col = "red")Both Models:
The QQ-plots for both models do not indicate that the residuals are normally distributed.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for both models.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$SOI_a_r_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Restricted Sociosexual Attitudes",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H3_M1, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M1 <- Boot(H3_M1, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M1) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -1.5075e-15 -0.0016338 0.057737 -0.0025904
TDDS_p_z 2.2087e-01 -0.0029972 0.055845 0.2187473
confint(boot_H3_M1) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1159650 0.1223375
TDDS_p_z 0.1132985 0.3316029
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, TDDS_p_z)We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and restricted sociosexual attitudes.
Based on the estimates, the standardized coefficient for pathogen disgust is around .22.
# Create a scatterplot with the regression line to visualize
plot(data$SOI_a_r_z, data$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Restricted Sociosexual Attitudes",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H3_M2, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
CRS_z <- as.numeric(data$CRS_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M2 <- Boot(H3_M2, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M2) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -2.4601e-16 0.0025578 0.053518 0.0036932
SOI_a_r_z 4.0934e-01 0.0016313 0.054043 0.4116882
confint(boot_H3_M2) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1133516 0.09850889
SOI_a_r_z 0.2938363 0.50917573
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z)We can see from the Scatter Plot that there is a positive trend in the relationship between restricted sociosexual attitudes and religiosity.
Based on the estimates, the standardized coefficient for restricted sociosexual attitudes is around .41.
The fact that pathogen disgust predicts a monogamous mating strategy and, in turn, a monogamous mating strategy predicts religiosity indicates that a monogamous mating strategy could potentially mediate the relationship between pathogen disgust and religiosity. Now we will run the mediation analysis.
# Fitting the full mediation model
H3_full <- lm(CRS_z ~ TDDS_p_z + SOI_a_r_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_full, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_full))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H3_full))
qqline(residuals(H3_full), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H3_full) TDDS_p_z SOI_a_r_z
1.051286 1.051286
# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_full <- Boot(H3_full, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_full) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -2.9133e-16 0.0017123 0.053076 0.0027606
TDDS_p_z 1.6040e-01 0.0021068 0.057003 0.1635359
SOI_a_r_z 3.7391e-01 0.0014924 0.056727 0.3760744
confint(boot_H3_full) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.10951763 0.09950991
TDDS_p_z 0.05084147 0.26724821
SOI_a_r_z 0.25700020 0.47854943
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H3_mediate_results <- mediate(model.m = H3_M1,
model.y = H3_full,
treat = "TDDS_p_z",
mediator = "SOI_a_r_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results of the mediation analysis
summary(H3_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.0826 0.0394 0.13 <2e-16 ***
ADE 0.1604 0.0515 0.26 0.008 **
Total Effect 0.2430 0.1304 0.35 <2e-16 ***
Prop. Mediated 0.3399 0.1682 0.67 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Mediation Analysis Results:
The ACME indicates that there is a significant indirect effect of pathogen disgust on religiosity through a monogamous mating strategy (B = .08, BCI = [.038, .13]).
The proportion of the total effect mediated (the ratio of the indirect effect to the direct effect) is .34. That is, 34% of the covariance between pathogen disgust and religiosity seems to be able to be accounted for by restricted sociosexual attitudes. Our inferential criterion for concluding full mediation was that this ratio would be at least .80, with the total effect greater than a B = .2. Even the upper bound of the bootstrapped confidence interval for this proportion mediated does not cross this threshold (BCI = [.164, .63]).
We will now re-test Hypothesis 1 while controlling for sex.
# Fitting the models
H1_M1_s <- lm(CONV_f_z ~ TDDS_p_z + sex, data = data)
H1_M2_s <- lm(CRS_z ~ CONV_f_z + sex, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_M1_s, which = 1) # For model 1plot(H1_M2_s, which = 1) # For model 2# Extract residuals and fitted values for model 1
residuals_m1 <- residuals(H1_M1_s)
fitted_m1 <- fitted(H1_M1_s)
# Plot residuals vs fitted values for model 1 with different colors for sex
plot(fitted_m1, residuals_m1, col = as.numeric(data$sex),
main = "Residuals vs Fitted Values (Model 1)",
xlab = "Fitted Values", ylab = "Residuals")
abline(h = 0, col = "red")
# Add a legend to the plot
legend("topright", legend = levels(data$sex), col = 1:length(levels(data$sex)), pch = 1)# Drop the residuals and fitted values from the environment
rm(residuals_m1, fitted_m1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_M1_s)) # For model 1hist(residuals(H1_M2_s)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H1_M1_s))
qqline(residuals(H1_M1_s), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H1_M2_s))
qqline(residuals(H1_M2_s), col = "red")Model 1:
Model 2:
We will conduct a bootstrapped tests of the predictors for both models.
Multicollinearity: Variance inflation factor (VIF) for both models
# Calculating the VIF for both models
vif(H1_M1_s)TDDS_p_z sex
1.038981 1.038981
vif(H1_M2_s)CONV_f_z sex
1.004388 1.004388
# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M1_s <- Boot(H1_M1_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_M1_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 0.068890 -0.00144952 0.087173 0.071407
TDDS_p_z 0.018766 0.00419217 0.066166 0.021675
sexFemale -0.139225 0.00018887 0.121467 -0.139117
confint(boot_H1_M1_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1101433 0.23420257
TDDS_p_z -0.1080678 0.14762599
sexFemale -0.3726060 0.08793543
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, TDDS_p_z, sex)Based on the estimate and BCI for pathogen disgust (B = .018, BCI = [-.114, .153]), controlling for sex does not change whether there is a positive relationship between pathogen disgust and traditionalism.
# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
CRS_z <- as.numeric(data$CRS_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M2_s <- Boot(H1_M2_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_M2_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -0.18812 0.0034793 0.075256 -0.18351
CONV_f_z 0.32702 0.0035275 0.051563 0.33122
sexFemale 0.38019 -0.0022518 0.109868 0.37625
confint(boot_H1_M2_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.3465750 -0.04560402
CONV_f_z 0.2115510 0.42600748
sexFemale 0.1722696 0.59940976
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z, sex)The fact that pathogen disgust does not predict traditionalism indicates that traditionalism could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H1_full_s <- lm(CRS_z ~ TDDS_p_z + CONV_f_z + sex, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_full_s, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_full_s))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H1_full_s))
qqline(residuals(H1_full_s), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a tendency for higher frequencies of residuals at just around 1 SD.
We will conduct bootstrapping for confidence intervals to make inferences.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H1_full_s)TDDS_p_z CONV_f_z sex
1.039334 1.004730 1.043866
# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_full_s <- Boot(H1_full_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_full_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -0.14722 0.00189747 0.074983 -0.14490
TDDS_p_z 0.21242 -0.00055277 0.057913 0.21240
CONV_f_z 0.32316 0.00279282 0.049039 0.32841
sexFemale 0.29753 -0.00086310 0.108512 0.29387
confint(boot_H1_full_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.29755901 -0.005976442
TDDS_p_z 0.09737364 0.322631503
CONV_f_z 0.21168286 0.409014822
sexFemale 0.09173913 0.515949990
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z, TDDS_p_z, sex)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H1_s_mediate_results <- mediate(model.m = H1_M1_s,
model.y = H1_full_s,
treat = "TDDS_p_z",
mediator = "CONV_f_z",
covariates = "sex",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results for the mediation analysis
summary(H1_s_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
(Inference Conditional on the Covariate Values Specified in `covariates')
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.00606 -0.03198 0.05 0.790
ADE 0.21242 0.09856 0.32 <2e-16 ***
Total Effect 0.21849 0.10509 0.34 0.002 **
Prop. Mediated 0.02776 -0.20174 0.25 0.788
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through traditionalism (B = .006, BCI = [-.032, .05]).
This result is inconsistent with both the weak and strong version version of Hypothesis 1. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by traditionalism, even when controlling for sex.
We will now re-test Hypothesis 2 while controlling for sex.
# Fitting the models
H2_M1_s <- lm(GENE_p_z ~ TDDS_p_z + sex, data = data)
H2_M2_s <- lm(CRS_z ~ GENE_p_z + sex, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_M1_s, which = 1) # For model 1plot(H2_M2_s, which = 1) # For model 2# Extract residuals and fitted values for model 1
residuals_m1 <- residuals(H2_M1_s)
fitted_m1 <- fitted(H2_M1_s)
# Plot residuals vs fitted values for model 1 with different colors for sex
plot(fitted_m1, residuals_m1, col = as.numeric(data$sex),
main = "Residuals vs Fitted Values (Model 1)",
xlab = "Fitted Values", ylab = "Residuals")
abline(h = 0, col = "red")
# Add a legend to the plot
legend("topright", legend = levels(data$sex), col = 1:length(levels(data$sex)), pch = 1)# Drop the residuals and fitted values from the environment
rm(residuals_m1, fitted_m1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_M1_s)) # For model 1hist(residuals(H2_M2_s)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H2_M1_s))
qqline(residuals(H2_M1_s), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H2_M2_s))
qqline(residuals(H2_M2_s), col = "red")Model 1:
Model 2:
We will conduct a bootstrapped tests of the predictors for both models.
Multicollinearity: Variance inflation factor (VIF) for both models
# Calculating the VIF for both models
vif(H2_M1_s)TDDS_p_z sex
1.038981 1.038981
vif(H2_M2_s)GENE_p_z sex
1.013159 1.013159
# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M1_s <- Boot(H2_M1_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_M1_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 0.119846 0.00029527 0.084531 0.117804
TDDS_p_z 0.037901 -0.00086675 0.060272 0.035864
sexFemale -0.242207 -0.00040661 0.117048 -0.242947
confint(boot_H2_M1_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.03909517 0.29625715
TDDS_p_z -0.08327969 0.15957373
sexFemale -0.46669146 -0.01264589
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, TDDS_p_z, sex)Based on the estimate and BCI for pathogen disgust (B = .038, BCI = [-.084, .157]), while controlling for sex there is not a relationship between pathogen disgust and out-group avoidance.
# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
CRS_z <- as.numeric(data$CRS_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2_s <- Boot(H2_M2_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_M2_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -0.17669 0.00383719 0.080908 -0.174657
GENE_p_z 0.08816 0.00058913 0.058323 0.089972
sexFemale 0.35710 -0.00291569 0.118240 0.350989
confint(boot_H2_M2_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.33163460 -0.01089098
GENE_p_z -0.03812488 0.19549662
sexFemale 0.13230476 0.58667840
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z, sex)The fact that pathogen disgust does not predict out-group avoidance and that out-group avoidance does not predict religiosity indicates that traditionalism could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H2_full_s <- lm(CRS_z ~ TDDS_p_z + GENE_p_z + sex, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_full_s, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_full_s))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H2_full_s))
qqline(residuals(H2_full_s), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a tendency for higher frequencies of residuals at just above -1 SD and just above 1 SD.
We will conduct bootstrapping for confidence intervals to make inferences.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H2_full_s)TDDS_p_z GENE_p_z sex
1.040438 1.014580 1.053911
# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_full_s <- Boot(H2_full_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_full_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -0.134569 0.00215572 0.080493 -0.131898
TDDS_p_z 0.215448 0.00086485 0.061324 0.214958
GENE_p_z 0.080198 0.00054884 0.056106 0.080718
sexFemale 0.271962 -0.00137026 0.117763 0.263523
confint(boot_H2_full_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.29354508 0.0244336
TDDS_p_z 0.09706498 0.3370680
GENE_p_z -0.03259225 0.1865872
sexFemale 0.05112016 0.5181156
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z, TDDS_p_z, sex)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_s_mediate_results <- mediate(model.m = H2_M1_s,
model.y = H2_full_s,
treat = "TDDS_p_z",
mediator = "GENE_p_z",
covariates = "sex",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results for the mediation analysis
summary(H2_s_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
(Inference Conditional on the Covariate Values Specified in `covariates')
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.00304 -0.00774 0.02 0.588
ADE 0.21545 0.10085 0.33 0.002 **
Total Effect 0.21849 0.10509 0.34 0.002 **
Prop. Mediated 0.01391 -0.04180 0.10 0.586
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through out-group avoidance (B = .003, BCI = [-.007, .02]).
This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by out-group avoidance, even when controlling for sex.
We will now re-test Hypothesis 3 while controlling for sex.
# Fitting the models
H3_M1_s <- lm(SOI_a_r_z ~ TDDS_p_z + sex, data = data)
H3_M2_s <- lm(CRS_z ~ SOI_a_r_z + sex, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_M1_s, which = 1) # For model 1plot(H3_M2_s, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_M1_s)) # For model 1hist(residuals(H3_M2_s)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H3_M1_s))
qqline(residuals(H3_M1_s), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H3_M2_s))
qqline(residuals(H3_M2_s), col = "red")Both Models:
We will conduct a bootstrapped tests of the predictors for both models.
Multicollinearity: Variance inflation factor (VIF) for both models
# Calculating the VIF for both models
vif(H3_M1_s)TDDS_p_z sex
1.038981 1.038981
vif(H3_M2_s)SOI_a_r_z sex
1.038623 1.038623
# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M1_s <- Boot(H3_M1_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M1_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -0.15403 -3.7669e-05 0.080601 -0.15429
TDDS_p_z 0.19067 -2.4728e-03 0.057487 0.18827
sexFemale 0.31129 -3.3652e-03 0.115848 0.30770
confint(boot_H3_M1_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.30713596 0.004964965
TDDS_p_z 0.07973518 0.305882636
sexFemale 0.09905334 0.546865791
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, TDDS_p_z, sex)# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
CRS_z <- as.numeric(data$CRS_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M2_s <- Boot(H3_M2_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M2_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -0.092212 0.0034832 0.078908 -0.08793
SOI_a_r_z 0.391337 0.0023672 0.054617 0.39386
sexFemale 0.186359 -0.0020876 0.114311 0.18463
confint(boot_H3_M2_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.24883828 0.06126489
SOI_a_r_z 0.27960214 0.49433536
sexFemale -0.02594773 0.40381904
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z, sex)# Fitting the full mediation model
H3_full_s <- lm(CRS_z ~ TDDS_p_z + SOI_a_r_z + sex, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_full_s, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_full_s))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H3_full_s))
qqline(residuals(H3_full_s), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not quite normally distributed.
We will conduct bootstrapping for confidence intervals to make inferences.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H3_full_s) TDDS_p_z SOI_a_r_z sex
1.078165 1.077794 1.065178
# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_full_s <- Boot(H3_full_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_full_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -0.069064 0.0022568 0.078290 -0.066787
TDDS_p_z 0.149296 0.0019682 0.058524 0.152164
SOI_a_r_z 0.362880 0.0020451 0.057094 0.365231
sexFemale 0.139577 -0.0013943 0.113992 0.137545
confint(boot_H3_full_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.22670927 0.07918258
TDDS_p_z 0.03465778 0.25852694
SOI_a_r_z 0.25101165 0.46914310
sexFemale -0.07659668 0.37316434
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z, TDDS_p_z, sex)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H3_s_mediate_results <- mediate(model.m = H3_M1_s,
model.y = H3_full_s,
treat = "TDDS_p_z",
mediator = "SOI_a_r_z",
covariates = "sex",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results for the mediation analysis
summary(H3_s_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
(Inference Conditional on the Covariate Values Specified in `covariates')
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.0692 0.0290 0.12 0.002 **
ADE 0.1493 0.0366 0.26 0.012 *
Total Effect 0.2185 0.1051 0.34 0.002 **
Prop. Mediated 0.3167 0.1444 0.69 0.004 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Mediation Analysis Results:
The ACME indicates that there is a significant indirect effect of pathogen disgust on religiosity through a monogamous mating strategy (B = .069, BCI = [.027, .12]).
The proportion of the total effect mediated (the ratio of the indirect effect to the direct effect) is .32. That is, 32% of the covariance between pathogen disgust and religiosity seems to be able to be accounted for by restricted sociosexual attitudes. Our inferential criterion for concluding full mediation was that this ratio would be at least .80, with the total effect greater than a B = .2. Even the upper bound of the bootstrapped confidence interval for this proportion mediated does not cross this threshold (BCI = [.136, .62]).
Now we will retest each hypothesis with religious importance as the outcome variable, as a check for robustness.
As mentioned above, to test Hypothesis 1, we must first establish that pathogen disgust affects adherence to traditional practices (i.e., CONV_f_z) and that adherence to traditional practices influences religious importance.
# Fitting the models
H1_M1_RI <- lm(CONV_f_z ~ TDDS_p_z, data = data)
H1_M2_RI <- lm(religious_importance_z ~ CONV_f_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_M1_RI, which = 1) # For model 1plot(H1_M2_RI, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_M1_RI)) # For model 1hist(residuals(H1_M2_RI)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H1_M1_RI))
qqline(residuals(H1_M1_RI), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H1_M2_RI))
qqline(residuals(H1_M2_RI), col = "red")Model 1:
Model 2:
The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$CONV_f_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Traditionalism",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H1_M1_RI, col = "black", lwd = 1)# Summarizing the model
summary(H1_M1_RI)
Call:
lm(formula = CONV_f_z ~ TDDS_p_z, data = data)
Residuals:
Min 1Q Median 3Q Max
-3.4517 -0.6429 0.0402 0.6287 2.8524
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.379e-17 5.893e-02 0.000 1.000
TDDS_p_z 5.260e-03 5.903e-02 0.089 0.929
Residual standard error: 1.002 on 287 degrees of freedom
Multiple R-squared: 2.766e-05, Adjusted R-squared: -0.003457
F-statistic: 0.00794 on 1 and 287 DF, p-value: 0.9291
We can see from the Scatter Plot that there is no relationship between pathogen disgust and traditionalism to speak of (B = .005, t(287) = .089, p = .929), just as with the correlation above.
# Create a scatterplot with the regression line to visualize
plot(data$CONV_f_z, data$religious_importance_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Traditionalism",
ylab = "Religious Importance",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H1_M2_RI, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M2_RI <- Boot(H1_M2_RI, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_M2_RI) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 1.7368e-16 0.0029611 0.055114 0.0011843
CONV_f_z 3.3199e-01 0.0044058 0.052992 0.3353958
confint(boot_H1_M2_RI) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1068043 0.1145093
CONV_f_z 0.2189193 0.4346032
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, religious_importance_z)We can see from the Scatter Plot that there is a positive relationship between traditionalism and religious importance.
Based on the estimates, the standardized coefficient for traditionalism is around .33.
The fact that pathogen disgust does not predict traditionalism indicates that traditionalism could not mediate the relationship between pathogen disgust and religious importance, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H1_full_RI <- lm(religious_importance_z ~ TDDS_p_z + CONV_f_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_full_RI, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_full_RI))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H1_full_RI))
qqline(residuals(H1_full_RI), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just below 1 SD.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H1_full_RI)TDDS_p_z CONV_f_z
1.000028 1.000028
# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_full_RI <- Boot(H1_full_RI, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_full_RI) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 1.8511e-16 0.00205054 0.053908 0.002309
TDDS_p_z 2.1682e-01 -0.00097327 0.057490 0.214199
CONV_f_z 3.3085e-01 0.00389128 0.050288 0.334003
confint(boot_H1_full_RI) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1082920 0.1108538
TDDS_p_z 0.1105920 0.3414131
CONV_f_z 0.2227669 0.4216086
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, religious_importance_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H1_RI_mediate_results <- mediate(model.m = H1_M1_RI,
model.y = H1_full_RI,
treat = "TDDS_p_z",
mediator = "CONV_f_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results for the mediation analysis
summary(H1_RI_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.00174 -0.03650 0.05 0.96
ADE 0.21682 0.09311 0.33 <2e-16 ***
Total Effect 0.21856 0.10172 0.34 <2e-16 ***
Prop. Mediated 0.00796 -0.22524 0.24 0.96
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs not containing zero, both pathogen disgust and traditionalism are significantly positively related to religious importance.
The coefficients for pathogen disgust (.22) and traditionalism (.33) are similar to the model with the CRS.
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religious importance through traditionalism (B = .002, BCI = [-.038, .04]).
This result is inconsistent with both the weak and strong version version of Hypothesis 1. That is, there is no evidence here that the relationship between pathogen disgust and religious importance is mediated by traditionalism.
As mentioned above, to test Hypothesis 2, we must first establish that pathogen disgust affects out-group avoidance (i.e., GENE_p_z) and that out-group avoidance influences religious importance.
# Fitting the models
H2_M1_RI <- lm(GENE_p_z ~ TDDS_p_z, data = data)
H2_M2_RI <- lm(religious_importance_z ~ GENE_p_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_M1_RI, which = 1) # For model 1plot(H2_M2_RI, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_M1_RI)) # For model 1hist(residuals(H2_M2_RI)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H2_M1_RI))
qqline(residuals(H2_M1_RI), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H2_M2_RI))
qqline(residuals(H2_M2_RI), col = "red")Model 1:
Model 2:
The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 or 1.5 SD.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$GENE_p_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Out-Group Avoidance (GENE_p)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M1_RI, col = "black", lwd = 1)# Summarizing the model
summary(H2_M1_RI)
Call:
lm(formula = GENE_p_z ~ TDDS_p_z, data = data)
Residuals:
Min 1Q Median 3Q Max
-1.9107 -0.8800 0.0804 0.5946 4.0239
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.976e-16 5.892e-02 0.000 1.000
TDDS_p_z 1.440e-02 5.902e-02 0.244 0.807
Residual standard error: 1.002 on 287 degrees of freedom
Multiple R-squared: 0.0002075, Adjusted R-squared: -0.003276
F-statistic: 0.05956 on 1 and 287 DF, p-value: 0.8074
We can see from the Scatter Plot that there is no relationship between pathogen disgust and out-group avoidance to speak of (B = .014, t(287) = .244, p = .807)—same as the model above.
# Create a scatterplot with the regression line to visualize
plot(data$GENE_p_z, data$religious_importance_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Out-Group Avoidance (GENE_p)",
ylab = "Religious Importance",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M2_RI, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2_RI <- Boot(H2_M2_RI, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_M2_RI) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 1.4098e-16 0.00271370 0.059242 0.0014941
GENE_p_z 7.4934e-02 -0.00037603 0.060076 0.0762719
confint(boot_H2_M2_RI) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.10931327 0.1249488
GENE_p_z -0.04800801 0.1879390
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, religious_importance_z)We can see from the Scatter Plot that there is a positive trend in the relationship between out-group avoidance and religious importance.
Based on the estimates, the standardized coefficient for traditionalism is around .07.
The fact that pathogen disgust does not predict out-group avoidance and, in turn, out-group avoidance does not predict religious importance indicates that out-group avoidance could not mediate the relationship between pathogen disgust and religious importance, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H2_full_RI <- lm(religious_importance_z ~ TDDS_p_z + GENE_p_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_full_RI, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_full_RI))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H2_full_RI))
qqline(residuals(H2_full_RI), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H2_full_RI)TDDS_p_z GENE_p_z
1.000208 1.000208
# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_full_RI <- Boot(H2_full_RI, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_full_RI) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 1.5313e-16 0.00183362 0.057942 0.0011976
TDDS_p_z 2.1752e-01 0.00041664 0.059876 0.2144363
GENE_p_z 7.1801e-02 -0.00026517 0.057177 0.0725622
confint(boot_H2_full_RI) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.10682993 0.1247495
TDDS_p_z 0.10889408 0.3438263
GENE_p_z -0.04473746 0.1872367
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, religious_importance_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_RI_mediate_results <- mediate(model.m = H2_M1_RI,
model.y = H2_full_RI,
treat = "TDDS_p_z",
mediator = "GENE_p_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results of the mediation analysis
summary(H2_RI_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.00103 -0.00929 0.02 0.84
ADE 0.21752 0.10126 0.34 <2e-16 ***
Total Effect 0.21856 0.10172 0.34 <2e-16 ***
Prop. Mediated 0.00473 -0.04825 0.08 0.84
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs, pathogen disgust, but not out-group avoidance, is significantly positively related to religious importance.
The coefficients for pathogen disgust (.22) and out-group avoidance (.07) are virtually identical in the model with CRS.
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religious importance through out-group avoidance (B = .001, BCI = [-.009, .01]).
This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religious importance is mediated by out-group avoidance.
As mentioned above, to test Hypothesis 3, we must first establish that pathogen disgust affects one’s mating strategy (i.e., make it more monogamous; as measured by more restricted sociosexuality; SOI_a_r_z) and that a monogamous mating strategy influences religious importance.
# Fitting the models
H3_M1_RI <- lm(SOI_a_r_z ~ TDDS_p_z, data = data)
H3_M2_RI <- lm(religious_importance_z ~ SOI_a_r_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_M1_RI, which = 1) # For model 1plot(H3_M2_RI, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_M1_RI)) # For model 1hist(residuals(H3_M2_RI)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H3_M1_RI))
qqline(residuals(H3_M1_RI), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H3_M2_RI))
qqline(residuals(H3_M2_RI), col = "red")Both Models:
The QQ-plots for both models do not indicate that the residuals are normally distributed.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for both models.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$SOI_a_r_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Restricted Sociosexual Attitudes",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H3_M1_RI, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M1_RI <- Boot(H3_M1_RI, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M1_RI) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -1.5075e-15 -0.0016338 0.057737 -0.0025904
TDDS_p_z 2.2087e-01 -0.0029972 0.055845 0.2187473
confint(boot_H3_M1_RI) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1159650 0.1223375
TDDS_p_z 0.1132985 0.3316029
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, TDDS_p_z)We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and restricted sociosexual attitudes.
Based on the estimates, the standardized coefficient for pathogen disgust is around .22.
# Create a scatterplot with the regression line to visualize
plot(data$SOI_a_r_z, data$religious_importance_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Restricted Sociosexual Attitudes",
ylab = "Religious Importance",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H3_M2_RI, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M2_RI <- Boot(H3_M2_RI, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M2_RI) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 7.4203e-16 0.0024378 0.054107 0.0011108
SOI_a_r_z 3.8592e-01 0.0025175 0.055091 0.3877059
confint(boot_H3_M2_RI) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1054531 0.1082780
SOI_a_r_z 0.2755024 0.4862312
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, religious_importance_z)We can see from the Scatter Plot that there is a positive trend in the relationship between restricted sociosexual attitudes and religious importance.
Based on the estimates, the standardized coefficient for restricted sociosexual attitudes is around .39.
The fact that pathogen disgust predicts a monogamous mating strategy and, in turn, a monogamous mating strategy predicts religious importance indicates that a monogamous mating strategy could potentially mediate the relationship between pathogen disgust and religious importance. Now we will run the mediation analysis.
# Fitting the full mediation model
H3_full_RI <- lm(religious_importance_z ~ TDDS_p_z + SOI_a_r_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_full_RI, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_full_RI))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H3_full_RI))
qqline(residuals(H3_full_RI), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H3_full_RI) TDDS_p_z SOI_a_r_z
1.051286 1.051286
# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_full_RI <- Boot(H3_full_RI, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_full_RI) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 7.0243e-16 0.0017189 0.053773 0.0011696
TDDS_p_z 1.4015e-01 0.0018125 0.057063 0.1415711
SOI_a_r_z 3.5496e-01 0.0021128 0.057777 0.3581158
confint(boot_H3_full_RI) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1063690 0.1042462
TDDS_p_z 0.0341517 0.2534335
SOI_a_r_z 0.2322735 0.4558705
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, religious_importance_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H3_RI_mediate_results <- mediate(model.m = H3_M1_RI,
model.y = H3_full_RI,
treat = "TDDS_p_z",
mediator = "SOI_a_r_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results of the mediation analysis
summary(H3_RI_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.0784 0.0371 0.13 <2e-16 ***
ADE 0.1402 0.0255 0.25 0.016 *
Total Effect 0.2186 0.1017 0.34 <2e-16 ***
Prop. Mediated 0.3587 0.1815 0.78 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Mediation Analysis Results:
The ACME indicates that there is a significant indirect effect of pathogen disgust on religious importance through a monogamous mating strategy (B = .078, BCI = [.037, .13]).
The proportion of the total effect mediated (the ratio of the indirect effect to the direct effect) is .36. That is, 36% of the covariance between pathogen disgust and religious importance seems to be able to be accounted for by restricted sociosexual attitudes. Our inferential criterion for concluding full mediation was that this ratio would be at least .80, with the total effect greater than a B = .2. Even the upper bound of the bootstrapped confidence interval for this proportion mediated does not cross this threshold (BCI = [.176, .71]).
In the Internal Reliability analysis of the Preferences Subscale of SFGENE-7, the drop alpha statistics indicated that removing Item 1 from the subscale would improve the internal reliability from an alpha of around .5 to around .8. Given that a lack of internal reliability may attenuate relationships between variables in correlational analyses, I will now conduct the same analysis for Hypothesis 2 but with Item 1 removed.
As mentioned above, to test Hypothesis 2, we must first establish that pathogen disgust affects out-group avoidance (i.e., GENE_p_2_z) and that out-group avoidance influences religiosity.
# Fitting the models
H2_M1_2 <- lm(GENE_p_2_z ~ TDDS_p_z, data = data)
H2_M2_2 <- lm(CRS_z ~ GENE_p_2_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_M1_2, which = 1) # For model 1plot(H2_M2_2, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_M1_2)) # For model 1hist(residuals(H2_M2_2)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H2_M1_2))
qqline(residuals(H2_M1_2), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H2_M2_2))
qqline(residuals(H2_M2_2), col = "red")Both Models:
The histograms and QQ-plots indicate that for both models the residuals are not normally distributed, although in different ways.
We will need to bootstrapping for testing the predictors for these models.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$GENE_p_2_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Out-Group Avoidance (GENE_p without Item 1)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M1_2, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
GENE_p_2_z <- as.numeric(data$GENE_p_2_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M1_2 <- Boot(H2_M1_2, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_M1_2) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -1.2029e-16 -0.00018739 0.059877 -0.0018552
TDDS_p_z -2.5360e-02 0.00068765 0.059533 -0.0261568
confint(boot_H2_M1_2) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1102701 0.12771218
TDDS_p_z -0.1342644 0.09246806
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_2_z, TDDS_p_z)We can see from the Scatter Plot that there is no relationship between pathogen disgust and out-group avoidance to speak of (B = -.025, BCI = [-.133, .106])—just as with the analysis above.
# Create a scatterplot with the regression line to visualize
plot(data$GENE_p_2_z, data$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Out-Group Avoidance (GENE_p without Item 1)",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M2_2, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
GENE_p_2_z <- as.numeric(data$GENE_p_2_z)
CRS_z <- as.numeric(data$CRS_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2_2 <- Boot(H2_M2_2, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_M2_2) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -8.6789e-16 2.2888e-03 0.059385 0.00088737
GENE_p_2_z 1.3388e-05 -2.0226e-05 0.062432 0.00152762
confint(boot_H2_M2_2) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1153198 0.1206342
GENE_p_2_z -0.1212795 0.1206176
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_2_z, CRS_z)We can see from the Scatter Plot that there is a positive trend in the relationship between out-group avoidance and religiosity.
Based on the estimates, the standardized coefficient for out-group avoidance is around B = .00001 (BCI = [-.119, .123]).
The fact that pathogen disgust does not predict out-group avoidance and, in turn, out-group avoidance does not predict religiosity indicates that out-group avoidance could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H2_2_full <- lm(CRS_z ~ TDDS_p_z + GENE_p_2_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_2_full, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_2_full))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H2_2_full))
qqline(residuals(H2_2_full), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H2_2_full) TDDS_p_z GENE_p_2_z
1.000644 1.000644
# Adding the variables to the broader environment so the Boot function will run
GENE_p_2_z <- as.numeric(data$GENE_p_2_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_2_full <- Boot(H2_2_full, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_2_full) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -8.5427e-16 0.00137003 0.057776 0.0028594
TDDS_p_z 2.4314e-01 0.00028211 0.059630 0.2416846
GENE_p_2_z 6.1794e-03 -0.00026584 0.057790 0.0055259
confint(boot_H2_2_full) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1188244 0.1086140
TDDS_p_z 0.1337811 0.3610994
GENE_p_2_z -0.1023601 0.1179598
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_2_z, CRS_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_2_mediate_results <- mediate(model.m = H2_M1_2,
model.y = H2_2_full,
treat = "TDDS_p_z",
mediator = "GENE_p_2_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results of the mediation analysis
summary(H2_2_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME -0.000157 -0.006777 0.01 0.92
ADE 0.243143 0.129699 0.35 <2e-16 ***
Total Effect 0.242987 0.130352 0.35 <2e-16 ***
Prop. Mediated -0.000645 -0.034402 0.05 0.92
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs, pathogen disgust, but not out-group avoidance, is significantly positively related to religiosity.
The coefficients for pathogen disgust (.24) and out-group avoidance (.006) are similar to their values with Item 1 included.
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through out-group avoidance (B = -.0002, BCI = [-.006, .01]).
This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by out-group avoidance.
Now we will rerun the analyses after ensuring that covariance between our mediating variables could not account for the results above. That is, we will regress each of the mediator on the other two and save standardized residuals. Then, we will run the mediation analyses again. This will help to explore whether any mediation effects (or null effects) in our previous analyses could have been due to covariation between our mediating variables. First, let us calculate our standardized residuals.
# Create the regression models for standardized residuals
# For traditionalism
res_model_CONV_f <- lm(CONV_f ~ SOI_a_r + GENE_p, data = data)
# For out-group avoidance
res_model_GENE_p <- lm(GENE_p ~ CONV_f + SOI_a_r, data = data)
# For sexual strategies
res_model_SOI_a_r <- lm(SOI_a_r ~ GENE_p + CONV_f, data = data)
# Add the standardized residuals to the new data frame, aligning them with the full rows
# Initializing new variables with NA values first
data$res_CONV_f <- NA # Initialize new variable with NAs for traditionalism
data$res_GENE_p <- NA # Initialize new variable with NAs for out-group avoidance
data$res_SOI_a_r <- NA
# Add the standardized residuals to the res_CONV_f variable
data$res_CONV_f <- rstandard(res_model_CONV_f)
# Add the standardized residuals to the res_GENE_p variable
data$res_GENE_p <- rstandard(res_model_GENE_p)
# Add the standardized residuals to the res_SOI_a_r variable
data$res_SOI_a_r <- rstandard(res_model_SOI_a_r)
# Ensure that the variables are plain numeric vectors
data$res_CONV_f <- as.numeric(data$res_CONV_f)
data$res_GENE_p <- as.numeric(data$res_GENE_p)
data$res_SOI_a_r <- as.numeric(data$res_SOI_a_r)As mentioned above, to test Hypothesis 1, we must first establish that pathogen disgust affects adherence to traditional practices (in this case the standardized residuals version; res_CONV_f) and that adherence to traditional practices influences religiosity.
# Fitting the models
H1_M1_res <- lm(res_CONV_f ~ TDDS_p_z, data = data)
H1_M2_res <- lm(CRS_z ~ res_CONV_f, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_M1_res, which = 1) # For model 1plot(H1_M2_res, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_M1_res)) # For model 1hist(residuals(H1_M2_res)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H1_M1_res))
qqline(residuals(H1_M1_res), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H1_M2_res))
qqline(residuals(H1_M2_res), col = "red")Model 1:
Model 2:
The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$res_CONV_f,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Traditionalism Standardized Residuals",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H1_M1_res, col = "black", lwd = 1)# Summarizing the model
summary(H1_M1_res)
Call:
lm(formula = res_CONV_f ~ TDDS_p_z, data = data)
Residuals:
Min 1Q Median 3Q Max
-3.5869 -0.6302 0.0261 0.6106 2.9054
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0002967 0.0590112 -0.005 0.996
TDDS_p_z -0.0272508 0.0591136 -0.461 0.645
Residual standard error: 1.003 on 287 degrees of freedom
Multiple R-squared: 0.0007399, Adjusted R-squared: -0.002742
F-statistic: 0.2125 on 1 and 287 DF, p-value: 0.6452
We can see from the Scatter Plot that there is no relationship between pathogen disgust and traditionalism to speak of (B = -.027, t(287) = .461, p = .645), just as with the analysis above.
# Create a scatterplot with the regression line to visualize
plot(data$res_CONV_f, data$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Traditionalism Standardized Residuals",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H1_M2_res, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
res_CONV_f <- as.numeric(data$res_CONV_f)
CRS_z <- as.numeric(data$CRS_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M2_res <- Boot(H1_M2_res, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_M2_res) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 7.5754e-05 0.0026181 0.056525 0.0028825
res_CONV_f 2.5535e-01 0.0034721 0.054933 0.2591966
confint(boot_H1_M2_res) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1156756 0.1095724
res_CONV_f 0.1345532 0.3565031
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(res_CONV_f, CRS_z)We can see from the Scatter Plot that there is a positive relationship between traditionalism and religiosity.
Based on the estimates, the standardized coefficient for traditionalism is around .25.
The fact that pathogen disgust does not predict traditionalism indicates that traditionalism could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H1_full_res <- lm(CRS_z ~ TDDS_p_z + res_CONV_f, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_full_res, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_full_res))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H1_full_res))
qqline(residuals(H1_full_res), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and around 1 SD.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H1_full_res) TDDS_p_z res_CONV_f
1.00074 1.00074
# Adding the variables to the broader environment so the Boot function will run
res_CONV_f <- as.numeric(data$res_CONV_f)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_full_res <- Boot(H1_full_res, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_full_res) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 7.7769e-05 0.00163568 0.054903 0.0027247
TDDS_p_z 2.5013e-01 -0.00066573 0.057291 0.2487006
res_CONV_f 2.6214e-01 0.00245588 0.050525 0.2645166
confint(boot_H1_full_res) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1114009 0.1040761
TDDS_p_z 0.1420620 0.3667127
res_CONV_f 0.1571168 0.3577343
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(res_CONV_f, CRS_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H1_mediate_results_res <- mediate(model.m = H1_M1_res,
model.y = H1_full_res,
treat = "TDDS_p_z",
mediator = "res_CONV_f",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results for the mediation analysis
summary(H1_mediate_results_res)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME -0.00714 -0.03622 0.03 0.63
ADE 0.25013 0.13001 0.36 <2e-16 ***
Total Effect 0.24299 0.13035 0.35 <2e-16 ***
Prop. Mediated -0.02940 -0.19132 0.12 0.63
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs not containing zero, both pathogen disgust and traditionalism are significantly positively related to religiosity.
The coefficients for pathogen disgust (.25) and traditionalism (.26) are close to how they are in the models above where they are stand-alone predictors. This suggests that the unique effects of pathogen disgust and traditionalism on religiosity over-and-above each other are independent effects.
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through traditionalism (B = -.007, BCI = [-.038, .03]).
This result is inconsistent with both the weak and strong version version of Hypothesis 1. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by traditionalism, even after removing covariation with the other potentially mediating variables.
As mentioned above, to test Hypothesis 2, we must first establish that pathogen disgust affects out-group avoidance (i.e., the standardized residual version of res_GENE_p) and that out-group avoidance influences religiosity.
# Fitting the models
H2_M1_res <- lm(res_GENE_p ~ TDDS_p_z, data = data)
H2_M2_res <- lm(CRS_z ~ res_GENE_p, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_M1_res, which = 1) # For model 1plot(H2_M2_res, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_M1_res)) # For model 1hist(residuals(H2_M2_res)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H2_M1_res))
qqline(residuals(H2_M1_res), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H2_M2_res))
qqline(residuals(H2_M2_res), col = "red")Model 1:
Model 2:
The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 or 1.5 SD.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$res_GENE_p,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Out-Group Avoidance (res_GENE_p) Standardized Residuals",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M1_res, col = "black", lwd = 1)# Summarizing the model
summary(H2_M1_res)
Call:
lm(formula = res_GENE_p ~ TDDS_p_z, data = data)
Residuals:
Min 1Q Median 3Q Max
-2.1221 -0.7581 -0.0587 0.7029 3.9597
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0001159 0.0590329 -0.002 0.998
TDDS_p_z -0.0175214 0.0591353 -0.296 0.767
Residual standard error: 1.004 on 287 degrees of freedom
Multiple R-squared: 0.0003058, Adjusted R-squared: -0.003177
F-statistic: 0.08779 on 1 and 287 DF, p-value: 0.7672
We can see from the Scatter Plot that there is no relationship between pathogen disgust and out-group avoidance to speak of (B = -.017, t(287) = -.296, p = .767), just as with the correlation above.
# Create a scatterplot with the regression line to visualize
plot(data$res_GENE_p, data$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Out-Group Avoidance (GENE_p)",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M2_res, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
res_GENE_p <- as.numeric(data$res_GENE_p)
CRS_z <- as.numeric(data$CRS_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2_res <- Boot(H2_M2_res, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_M2_res) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -1.1983e-06 2.4874e-03 0.059214 0.0011912
res_GENE_p -1.0335e-02 -9.9515e-06 0.061660 -0.0118285
confint(boot_H2_M2_res) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1148799 0.1187364
res_GENE_p -0.1243771 0.1173597
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(res_GENE_p, CRS_z)We can see from the Scatter Plot that there no relationship.
Based on the estimates, the standardized coefficient for out-group avoidance is around .01.
The fact that pathogen disgust does not predict out-group avoidance and, in turn, out-group avoidance does not predict religiosity indicates that out-group avoidance could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H2_full_res <- lm(CRS_z ~ TDDS_p_z + res_GENE_p, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_full_res, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_full_res))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H2_full_res))
qqline(residuals(H2_full_res), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H2_full_res) TDDS_p_z res_GENE_p
1.000306 1.000306
# Adding the variables to the broader environment so the Boot function will run
res_GENE_p <- as.numeric(data$res_GENE_p)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_full_res <- Boot(H2_full_res, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_full_res) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -7.0685e-07 0.00151630 0.057644 0.0022743
TDDS_p_z 2.4288e-01 0.00069823 0.059676 0.2411531
res_GENE_p -6.0965e-03 -0.00012461 0.058710 -0.0056310
confint(boot_H2_full_res) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1181186 0.1086985
TDDS_p_z 0.1325909 0.3610056
res_GENE_p -0.1214375 0.1088764
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(res_GENE_p, CRS_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_mediate_results_res <- mediate(model.m = H2_M1_res,
model.y = H2_full_res,
treat = "TDDS_p_z",
mediator = "res_GENE_p",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results of the mediation analysis
summary(H2_mediate_results_res)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.000107 -0.006742 0.01 0.87
ADE 0.242880 0.130981 0.36 <2e-16 ***
Total Effect 0.242987 0.130352 0.35 <2e-16 ***
Prop. Mediated 0.000440 -0.027246 0.04 0.87
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs, pathogen disgust, but not out-group avoidance, is significantly positively related to religiosity.
The coefficients for pathogen disgust (.24) and out-group avoidance (.-.006) are very similar in this model to the models above where they are stand-alone predictors.
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through out-group avoidance (B = .0001, BCI = [-.006, .01]).
This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by out-group avoidance, even when using the standardized residuals.
As mentioned above, to test Hypothesis 3, we must first establish that pathogen disgust affects one’s mating strategy (i.e., make it more monogamous; as measured by more restricted sociosexuality; in this case the standardized residuals version res_SOI_a_r) and that a monogamous mating strategy influences religiosity.
# Fitting the models
H3_M1_res <- lm(res_SOI_a_r ~ TDDS_p_z, data = data)
H3_M2_res <- lm(CRS_z ~ res_SOI_a_r, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_M1_res, which = 1) # For model 1plot(H3_M2_res, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_M1_res)) # For model 1hist(residuals(H3_M2_res)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H3_M1_res))
qqline(residuals(H3_M1_res), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H3_M2_res))
qqline(residuals(H3_M2_res), col = "red")Both Models:
The QQ-plots for both models do not indicate that the residuals are normally distributed.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for both models.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$res_SOI_a_r,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Restricted Sociosexual Attitudes",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H3_M1_res, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
res_SOI_a_r <- as.numeric(data$res_SOI_a_r)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M1_res <- Boot(H3_M1_res, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M1_res) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 4.5205e-05 -0.0015333 0.057547 -0.0027783
TDDS_p_z 2.2359e-01 -0.0035998 0.054417 0.2211374
confint(boot_H3_M1_res) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1134134 0.1188434
TDDS_p_z 0.1248918 0.3356576
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(res_SOI_a_r, TDDS_p_z)We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and restricted sociosexual attitudes.
Based on the estimates, the standardized coefficient for pathogen disgust is around .22.
# Create a scatterplot with the regression line to visualize
plot(data$res_SOI_a_r, data$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Restricted Sociosexual Attitudes",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H3_M2_res, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
res_SOI_a_r <- as.numeric(data$res_SOI_a_r)
CRS_z <- as.numeric(data$CRS_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M2_res <- Boot(H3_M2_res, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M2_res) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -1.6429e-05 0.0024569 0.054711 0.0038558
res_SOI_a_r 3.6343e-01 0.0020824 0.056629 0.3659118
confint(boot_H3_M2_res) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1119157 0.1047078
res_SOI_a_r 0.2401447 0.4617742
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(res_SOI_a_r, CRS_z)We can see from the Scatter Plot that there is a positive trend in the relationship between restricted sociosexual attitudes and religiosity.
Based on the estimates, the standardized coefficient for restricted sociosexual attitudes is around .36.
The fact that pathogen disgust predicts a monogamous mating strategy and, in turn, a monogamous mating strategy predicts religiosity indicates that a monogamous mating strategy could potentially mediate the relationship between pathogen disgust and religiosity. Now we will run the mediation analysis.
# Fitting the full mediation model
H3_full_res <- lm(CRS_z ~ TDDS_p_z + res_SOI_a_r, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_full_res, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_full_res))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H3_full_res))
qqline(residuals(H3_full_res), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H3_full_res) TDDS_p_z res_SOI_a_r
1.052405 1.052405
# Adding the variables to the broader environment so the Boot function will run
res_SOI_a_r <- as.numeric(data$res_SOI_a_r)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_full_res <- Boot(H3_full_res, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_full_res) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -1.4715e-05 0.0015772 0.054194 0.0017929
TDDS_p_z 1.7020e-01 0.0020435 0.058640 0.1728617
res_SOI_a_r 3.2552e-01 0.0020325 0.059172 0.3280743
confint(boot_H3_full_res) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1088316 0.1041097
TDDS_p_z 0.0562329 0.2816977
res_SOI_a_r 0.2050904 0.4359742
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(res_SOI_a_r, CRS_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H3_mediate_results_res <- mediate(model.m = H3_M1_res,
model.y = H3_full_res,
treat = "TDDS_p_z",
mediator = "res_SOI_a_r",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results of the mediation analysis
summary(H3_mediate_results_res)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.0728 0.0346 0.12 <2e-16 ***
ADE 0.1702 0.0603 0.27 0.006 **
Total Effect 0.2430 0.1304 0.35 <2e-16 ***
Prop. Mediated 0.2995 0.1471 0.60 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Mediation Analysis Results:
The ACME indicates that there is a significant indirect effect of pathogen disgust on religiosity through a monogamous mating strategy (B = .07, BCI = [.03, .12]).
The proportion of the total effect mediated (the ratio of the indirect effect to the direct effect) is .30. That is, 30% of the covariance between pathogen disgust and religiosity seems to be able to be accounted for by restricted sociosexual attitudes. Our inferential criterion for concluding full mediation was that this ratio would be at least .80, with the total effect greater than a B = .2. Even the upper bound of the bootstrapped confidence interval for this proportion mediated does not cross this threshold (BCI = [.13, .57]).
Ultimately, the results of reconducting the analysis with standardized residuals reinforces the results of the analysis thus far.
After receiving helpful suggestions from the reviewers at Religion, Brain and Behavior, we have decided to conduct a number of exploratory analyses that will disambiguate the findings here.
Before testing our hypotheses with only religious people, we will filter only those that indicated that they are religious
# Filter data to include only religious participants (excludes "None", "Anti-religious", "Athiest")
data_rel <- data %>%
filter(religious_affiliation %in% c("Christian", "Muslim", "Buddhist",
"Hindu", "Pagan", "Spiritual", "Spiritualist", "Deist"))
# Reversing the scoring for the SOI_a (saving it as SOI_a_r) and standardizing all variables
data_rel <- data_rel %>%
mutate(
SOI_a_r = (max(data_rel$SOI_a) + 1) - SOI_a,
SOI_a_r_z = scale(SOI_a_r),
CRS_z = scale(CRS),
religious_importance_z = scale(religious_importance),
GENE_p_z = scale(GENE_p),
GENE_p_2_z = scale(GENE_p_2),
CONV_f_z = scale(CONV_f),
TDDS_p_z = scale(TDDS_p)
)
# Ensure that the variables are plain numeric vectors
data_rel$TDDS_p_z <- as.numeric(data_rel$TDDS_p_z)
data_rel$CONV_f_z <- as.numeric(data_rel$CONV_f_z)
data_rel$CRS_z <- as.numeric(data_rel$CRS_z)
data_rel$SOI_a_r_z <- as.numeric(data_rel$SOI_a_r_z)
data_rel$religious_importance_z <- as.numeric(data_rel$religious_importance_z)
data_rel$GENE_p_z <- as.numeric(data_rel$GENE_p_z)
data_rel$GENE_p_2_z <- as.numeric(data_rel$GENE_p_2_z)
# Because the boot function was throwing a fit about this earlier
CRS_z <- as.numeric(data_rel$CRS_z)
TDDS_p_z <- as.numeric(data_rel$TDDS_p_z)Now, there is n = 141 participants for analysis here. All are religious.
Before moving on, we will first test the baseline model where we regress the dependent variable (CRS) on the independent variable (TDDS_p), because each mediation analysis requires that there is a significant positive relationship between these variables. We will also test each of the assumptions.
# Fitting the model
baseline_model_rel <- lm(CRS_z ~ TDDS_p_z, data = data_rel)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(baseline_model_rel, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(baseline_model_rel))# QQ-plot for normality of residuals
qqnorm(residuals(baseline_model_rel))
qqline(residuals(baseline_model_rel), col = "red")The histogram of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a tendency for higher frequencies of residuals at both just above -1 SD and just below 1 SD.
Because our inferential criteria do not require parametric tests, we should not be in too much trouble with a non-normality of residuals. We were relying on a parametric t-test to assess whether pathogen disgust is related to religiosity in the first place, however. Because we have failed the assumption of normality of residuals, we will use a bootstrapped test to make our inference here as well.
Multicollinearity: Not applicable becuase only one predictor
# Create a scatterplot with the regression line to visualize
plot(data_rel$TDDS_p_z, data_rel$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(baseline_model_rel, col = "black", lwd = 1)# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_baseline_model_rel <- Boot(baseline_model_rel, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_baseline_model_rel) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 2.3086e-16 -7.2771e-05 0.079804 -0.0022716
TDDS_p_z 3.1293e-01 1.9212e-03 0.086146 0.3212187
confint(boot_baseline_model_rel) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1471255 0.1619620
TDDS_p_z 0.1143230 0.4582468
# Removing the CRS_z and TDDS_p_z vectors in the environment that the Boot package required outside of the data.frame
rm(CRS_z, TDDS_p_z)We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and religiosity.
Based on the estimates, the standardized coefficient for pathogen disgust is around .31. This is slightly higher than for the whole sample, including non-religious people (with was .24).
This shows that the IV effects the DV, as in the second regression equation for testing mediation in Baron & Kenny (1986). To establish mediation, we will still need to establish that the IV affects the mediator and the mediator affects the DV. Before we test the indirect effect and proportion of variance explained by each effect for each mediation model, we will run these two models.
As mentioned above, to test Hypothesis 1, we must first establish that pathogen disgust affects adherence to traditional practices (i.e., CONV_f_z) and that adherence to traditional practices influences religiosity.
# Fitting the models
H1_M1_rel <- lm(CONV_f_z ~ TDDS_p_z, data = data_rel)
H1_M2_rel <- lm(CRS_z ~ CONV_f_z, data = data_rel)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_M1_rel, which = 1) # For model 1plot(H1_M2_rel, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_M1_rel)) # For model 1hist(residuals(H1_M2_rel)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H1_M1_rel))
qqline(residuals(H1_M1_rel), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H1_M2_rel))
qqline(residuals(H1_M2_rel), col = "red")Model 1:
Model 2:
The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just below 1 SD.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data_rel$TDDS_p_z, data_rel$CONV_f_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Traditionalism",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H1_M1_rel, col = "black", lwd = 1)# Summarizing the model
summary(H1_M1_rel)
Call:
lm(formula = CONV_f_z ~ TDDS_p_z, data = data_rel)
Residuals:
Min 1Q Median 3Q Max
-3.4199 -0.5498 -0.0864 0.6819 2.9068
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.454e-16 8.428e-02 0.000 1.000
TDDS_p_z -7.504e-02 8.458e-02 -0.887 0.377
Residual standard error: 1.001 on 139 degrees of freedom
Multiple R-squared: 0.005631, Adjusted R-squared: -0.001523
F-statistic: 0.7871 on 1 and 139 DF, p-value: 0.3765
We can see from the Scatter Plot that there is a slight negative relationship between pathogen disgust and traditionalism, but this relationship is not significant (B = -.07, t(139) = -.887, p = .377).
This negative result (while there is a null relationship in the population) could be an artifact of collider bias, where to be religious you could either be traditional or disgust sensitive and therefore in the population of religious people there is a slight negative relationship between these.
Although this precludes the possibility that adherence to tradition mediates the relationship between disgust and religiosity, we will continue with the analysis for completions sake.
# Create a scatterplot with the regression line to visualize
plot(data_rel$CONV_f_z, data_rel$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Traditionalism",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H1_M2_rel, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data_rel$CONV_f_z)
CRS_z <- as.numeric(data_rel$CRS_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M2_rel <- Boot(H1_M2_rel, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_M2_rel) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 2.6806e-16 -0.00086015 0.083509 -0.0012939
CONV_f_z 8.6179e-02 0.00724783 0.089536 0.0925694
confint(boot_H1_M2_rel) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1622856 0.1598681
CONV_f_z -0.1106986 0.2442643
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z)We can see from the Scatter Plot that there is a slight positive relationship between traditionalism and religiosity.
Based on the estimates, the standardized coefficient for traditionalism is around .09.
The fact that pathogen disgust does not predict traditionalism indicates that traditionalism could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H1_full_rel <- lm(CRS_z ~ TDDS_p_z + CONV_f_z, data = data_rel)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_full_rel, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_full_rel))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H1_full_rel))
qqline(residuals(H1_full_rel), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a tendency for higher frequencies of residuals at both just above -1 SD and just below 1 SD.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H1_full_rel)TDDS_p_z CONV_f_z
1.005663 1.005663
# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data_rel$CONV_f_z)
CRS_z <- as.numeric(data_rel$CRS_z)
TDDS_p_z <- as.numeric(data_rel$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_full_rel <- Boot(H1_full_rel, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_full_rel) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 2.4689e-16 -9.7424e-05 0.079554 -0.001796
TDDS_p_z 3.2121e-01 4.5498e-04 0.088147 0.325711
CONV_f_z 1.1028e-01 6.8818e-03 0.078386 0.119082
confint(boot_H1_full_rel) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.14873819 0.1608342
TDDS_p_z 0.11979446 0.4743443
CONV_f_z -0.06180088 0.2456673
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H1_mediate_results_rel <- mediate(model.m = H1_M1_rel,
model.y = H1_full_rel,
treat = "TDDS_p_z",
mediator = "CONV_f_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results for the mediation analysis
summary(H1_mediate_results_rel)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME -0.00828 -0.03816 0.02 0.550
ADE 0.32121 0.15388 0.48 0.002 **
Total Effect 0.31293 0.14954 0.47 0.002 **
Prop. Mediated -0.02644 -0.15190 0.08 0.548
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 141
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs not containing zero, pathogen disgust, but not traditionalism, is significantly positively related to religiosity.
The coefficients: for pathogen disgust (.32) and traditionalism (.11).
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through traditionalism (B = -.008, BCI = [-.044, .02]).
This result is inconsistent with both the weak and strong version version of Hypothesis 1. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by traditionalism, including when looking at religious people only.
As mentioned above, to test Hypothesis 2, we must first establish that pathogen disgust affects out-group avoidance (i.e., GENE_p_z) and that out-group avoidance influences religiosity.
# Fitting the models
H2_M1_rel <- lm(GENE_p_z ~ TDDS_p_z, data = data_rel)
H2_M2_rel <- lm(CRS_z ~ GENE_p_z, data = data_rel)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_M1_rel, which = 1) # For model 1plot(H2_M2_rel, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_M1_rel)) # For model 1hist(residuals(H2_M2_rel)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H2_M1_rel))
qqline(residuals(H2_M1_rel), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H2_M2_rel))
qqline(residuals(H2_M2_rel), col = "red")Model 1:
Model 2:
The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just below th mean and just below 1 SD.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data_rel$TDDS_p_z, data_rel$GENE_p_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Out-Group Avoidance (GENE_p)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M1_rel, col = "black", lwd = 1)# Summarizing the model
summary(H2_M1_rel)
Call:
lm(formula = GENE_p_z ~ TDDS_p_z, data = data_rel)
Residuals:
Min 1Q Median 3Q Max
-2.0571 -0.6126 -0.0489 0.6251 3.8059
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.021e-16 8.438e-02 0.000 1.000
TDDS_p_z 5.742e-02 8.468e-02 0.678 0.499
Residual standard error: 1.002 on 139 degrees of freedom
Multiple R-squared: 0.003297, Adjusted R-squared: -0.003873
F-statistic: 0.4598 on 1 and 139 DF, p-value: 0.4988
We can see from the Scatter Plot that there is no relationship between pathogen disgust and out-group avoidance to speak of (B = .057, t(139) = .678, p = .499).
# Create a scatterplot with the regression line to visualize
plot(data_rel$GENE_p_z, data_rel$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Out-Group Avoidance (GENE_p)",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M2_rel, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data_rel$GENE_p_z)
CRS_z <- as.numeric(data_rel$CRS_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2_rel <- Boot(H2_M2_rel, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_M2_rel) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 2.6622e-16 -0.00024138 0.084047 0.00042404
GENE_p_z -5.4154e-02 -0.00716174 0.088882 -0.05820085
confint(boot_H2_M2_rel) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1666378 0.1610149
GENE_p_z -0.2223690 0.1094601
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z)We can see from the Scatter Plot that there is a negative trend in the relationship between out-group avoidance and religiosity.
Based on the estimates, the standardized coefficient for out-group avoidance is around -.05.
The fact that pathogen disgust does not predict out-group avoidance and, in turn, out-group avoidance does not predict religiosity indicates that out-group avoidance could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H2_full_rel <- lm(CRS_z ~ TDDS_p_z + GENE_p_z, data = data_rel)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_full_rel, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_full_rel))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H2_full_rel))
qqline(residuals(H2_full_rel), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a tendency for higher frequencies of residuals at both just above -1 SD and just below 1 SD.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H2_full_rel)TDDS_p_z GENE_p_z
1.003308 1.003308
# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data_rel$GENE_p_z)
CRS_z <- as.numeric(data_rel$CRS_z)
TDDS_p_z <- as.numeric(data_rel$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_full_rel <- Boot(H2_full_rel, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_full_rel) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 2.4549e-16 -7.8376e-05 0.080112 -0.0015804
TDDS_p_z 3.1709e-01 1.0441e-04 0.087176 0.3234525
GENE_p_z -7.2362e-02 -4.9109e-03 0.079234 -0.0740654
confint(boot_H2_full_rel) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1502413 0.15947547
TDDS_p_z 0.1276348 0.46989918
GENE_p_z -0.2260344 0.07345632
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_mediate_results_rel <- mediate(model.m = H2_M1_rel,
model.y = H2_full_rel,
treat = "TDDS_p_z",
mediator = "GENE_p_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results of the mediation analysis
summary(H2_mediate_results_rel)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME -0.00416 -0.02253 0.02 0.708
ADE 0.31709 0.14960 0.47 0.002 **
Total Effect 0.31293 0.14954 0.47 0.002 **
Prop. Mediated -0.01328 -0.08283 0.06 0.710
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 141
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs, pathogen disgust, but not out-group avoidance, is significantly positively related to religiosity.
The coefficients: for pathogen disgust (.32) and out-group avoidance (-.07).
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through out-group avoidance (B = -.004, BCI = [-.023, .02]).
This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by out-group avoidance, even if only assessing religious participants.
As mentioned above, to test Hypothesis 3, we must first establish that pathogen disgust affects one’s mating strategy (i.e., make it more monogamous; as measured by more restricted sociosexuality; SOI_a_r_z) and that a monogamous mating strategy influences religiosity.
# Fitting the models
H3_M1_rel <- lm(SOI_a_r_z ~ TDDS_p_z, data = data_rel)
H3_M2_rel <- lm(CRS_z ~ SOI_a_r_z, data = data_rel)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_M1_rel, which = 1) # For model 1plot(H3_M2_rel, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_M1_rel)) # For model 1hist(residuals(H3_M2_rel)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H3_M1_rel))
qqline(residuals(H3_M1_rel), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H3_M2_rel))
qqline(residuals(H3_M2_rel), col = "red")Both Models:
The QQ-plots for both models do not indicate that the residuals are normally distributed.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for both models.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data_rel$TDDS_p_z, data_rel$SOI_a_r_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Restricted Sociosexual Attitudes",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H3_M1_rel, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data_rel$SOI_a_r_z)
TDDS_p_z <- as.numeric(data_rel$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M1_rel <- Boot(H3_M1_rel, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M1_rel) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 1.3407e-15 0.00077784 0.081682 0.0028432
TDDS_p_z 1.9972e-01 0.00086515 0.080392 0.1985761
confint(boot_H3_M1_rel) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1664251 0.1462323
TDDS_p_z 0.0464425 0.3473499
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, TDDS_p_z)We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and restricted sociosexual attitudes.
Based on the estimates, the standardized coefficient for pathogen disgust is around .19.
# Create a scatterplot with the regression line to visualize
plot(data_rel$SOI_a_r_z, data_rel$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Restricted Sociosexual Attitudes",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H3_M2_rel, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data_rel$SOI_a_r_z)
CRS_z <- as.numeric(data_rel$CRS_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M2_rel <- Boot(H3_M2_rel, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M2_rel) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -2.6447e-16 -0.0013189 0.078758 -0.0016614
SOI_a_r_z 3.8309e-01 -0.0043342 0.080185 0.3809897
confint(boot_H3_M2_rel) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1610559 0.1580375
SOI_a_r_z 0.2218491 0.5358196
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z)We can see from the Scatter Plot that there is a positive trend in the relationship between restricted sociosexual attitudes and religiosity.
Based on the estimates, the standardized coefficient for restricted sociosexual attitudes is around .38.
The fact that pathogen disgust predicts a monogamous mating strategy and, in turn, a monogamous mating strategy predicts religiosity indicates that a monogamous mating strategy could potentially mediate the relationship between pathogen disgust and religiosity. Now we will run the mediation analysis.
# Fitting the full mediation model
H3_full_rel <- lm(CRS_z ~ TDDS_p_z + SOI_a_r_z, data = data_rel)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_full_rel, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_full_rel))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H3_full_rel))
qqline(residuals(H3_full_rel), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a tendency for higher frequencies of residuals at both just below the mean and just below 1 SD.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H3_full_rel) TDDS_p_z SOI_a_r_z
1.041546 1.041546
# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data_rel$SOI_a_r_z)
CRS_z <- as.numeric(data_rel$CRS_z)
TDDS_p_z <- as.numeric(data_rel$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_full_rel <- Boot(H3_full_rel, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_full_rel) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -2.1680e-16 -0.0012018 0.076513 -0.0015274
TDDS_p_z 2.4625e-01 0.0023376 0.086831 0.2520226
SOI_a_r_z 3.3391e-01 -0.0043978 0.082107 0.3283799
confint(boot_H3_full_rel) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.14622859 0.1585309
TDDS_p_z 0.05404243 0.4148393
SOI_a_r_z 0.17970551 0.4927243
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H3_mediate_results_rel <- mediate(model.m = H3_M1_rel,
model.y = H3_full_rel,
treat = "TDDS_p_z",
mediator = "SOI_a_r_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results of the mediation analysis
summary(H3_mediate_results_rel)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.0667 0.0141 0.14 0.012 *
ADE 0.2462 0.0806 0.41 0.006 **
Total Effect 0.3129 0.1495 0.47 0.002 **
Prop. Mediated 0.2131 0.0479 0.53 0.014 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 141
Simulations: 1000
Multiple Regression Results:
Mediation Analysis Results:
The ACME indicates that there is a significant indirect effect of pathogen disgust on religiosity through a monogamous mating strategy (B = .067, BCI = [.013, .13]).
The proportion of the total effect mediated (the ratio of the indirect effect to the direct effect) is .21. That is, 21% of the covariance between pathogen disgust and religiosity seems to be able to be accounted for by restricted sociosexual attitudes. Our inferential criterion for concluding full mediation was that this ratio would be at least .80, with the total effect greater than a B = .2. Even the upper bound of the bootstrapped confidence interval for this proportion mediated does not cross this threshold (BCI = [.044, .54]).
Ultimately, the results of including only religious individuals in the analysis support the results thus far. They are in line with the weak version of the sexual strategies hypothesis only. In addition, the proportion mediated is slightly smaller when analyzing only religious individuals (21% vs. 34%). This suggests that while the sexual strategies hypothesis may explain some differences in religiosity between religious people, it is better at explaining individual differences in religiosity including those that are not explicitly religious. This may simply be a feature of restricted range, or it may suggest that the sexual strategies hypothesis can explain both why religious people differ in their “devotion” and why some people are religious and others are not.
BIS reactivity may encourage traditional cultural values (like religion) as a means of mitigating social behaviors (like sex) that are more likely to result in the spread of infectious disease. We have focused thus far on religion in this analysis as the cultural value influenced by BIS reactivity and the specific mechanisms through which this may be the case. We will now analyze whether sociosexual attitudes mediate the relation between sexual disgust and traditionalism.
We will first test the baseline model where we regress the dependent variable (CONV_f_z) on the independent variable (TDDS_s_z), because each mediation analysis requires that there is a significant positive relationship between these variables. We will also test each of the assumptions. We will first need to create the TDDS_s_z variable by standardizing TDDS_s.
# Creating the standardized sexual disgust variable
data$TDDS_s_z <- scale(data$TDDS_s)
# Ensure it is numeric so the boot function will run
data$TDDS_s_z <- as.numeric(data$TDDS_s_z)
# Fitting the model
baseline_model_SST <- lm(CONV_f_z ~ TDDS_s_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(baseline_model_SST, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(baseline_model_SST))# QQ-plot for normality of residuals
qqnorm(residuals(baseline_model_SST))
qqline(residuals(baseline_model_SST), col = "red")Multicollinearity: Not applicable becuase only one predictor
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_s_z, data$CONV_f_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Sexual Disgust",
ylab = "Traditionalism",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(baseline_model_SST, col = "black", lwd = 1)# Summarizing the model
summary(baseline_model_SST)
Call:
lm(formula = CONV_f_z ~ TDDS_s_z, data = data)
Residuals:
Min 1Q Median 3Q Max
-3.3762 -0.5975 0.0059 0.6273 2.7846
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.880e-17 5.842e-02 0.000 1.0000
TDDS_s_z 1.307e-01 5.852e-02 2.234 0.0263 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9931 on 287 degrees of freedom
Multiple R-squared: 0.01709, Adjusted R-squared: 0.01366
F-statistic: 4.989 on 1 and 287 DF, p-value: 0.02628
We can see from the Scatter Plot that there is a slight positive relationship between sexual disgust and religiosity.
Based on the estimates, the (significant) standardized coefficient for sexual disgust is around .131 (p = .026).
This shows that the IV effects the DV, as in the second regression equation for testing mediation in Baron & Kenny (1986). To establish mediation, we will still need to establish that the IV affects the mediator and the mediator affects the DV. Before we test the indirect effect and proportion of variance explained by each effect for each mediation model, we will run these two models.
As mentioned above, to test our mediaiton model, we must first establish that sexual disgust affects sociosexual attitudes (i.e., SOI_a_r_z) and that sociosexual attitudes influences traditionalism.
# Fitting the models
SST_M1 <- lm(SOI_a_r_z ~ TDDS_s_z, data = data)
SST_M2 <- lm(CONV_f_z ~ SOI_a_r_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(SST_M1, which = 1) # For model 1plot(SST_M2, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(SST_M1)) # For model 1hist(residuals(SST_M2)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(SST_M1))
qqline(residuals(SST_M1), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(SST_M2))
qqline(residuals(SST_M2), col = "red")Both Models:
We will conduct bootstrapping for Model 1 because of the ambivalent residuals plot above.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_s_z, data$SOI_a_r_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Restricted Sociosexual Attitudes",
ylab = "Sexual Disgust",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(SST_M1, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
TDDS_s_z <- as.numeric(data$TDDS_s_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_SST_M1 <- Boot(SST_M1, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_SST_M1) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -1.5856e-15 -0.0020961 0.046731 -0.001715
TDDS_s_z 5.8840e-01 0.0002019 0.042232 0.588610
confint(boot_SST_M1) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.09047559 0.09164489
TDDS_s_z 0.50214194 0.66703586
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, TDDS_s_z)# Create a scatterplot with the regression line to visualize
plot(data$SOI_a_r_z, data$CONV_f_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Restricted Sociosexual Attitudes",
ylab = "Traditionalism",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(SST_M2, col = "black", lwd = 1)# Summarizing the model
summary(SST_M2)
Call:
lm(formula = CONV_f_z ~ SOI_a_r_z, data = data)
Residuals:
Min 1Q Median 3Q Max
-3.5045 -0.6536 0.0191 0.6226 2.7662
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.739e-16 5.826e-02 0.000 1.0000
SOI_a_r_z 1.500e-01 5.836e-02 2.571 0.0107 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9904 on 287 degrees of freedom
Multiple R-squared: 0.02251, Adjusted R-squared: 0.0191
F-statistic: 6.609 on 1 and 287 DF, p-value: 0.01065
We can see from the Scatter Plot that there is a slight positive relationship between sociosexual attitudes and traditionalism.
# Fitting the full mediation model
SST_full <- lm(CONV_f_z ~ TDDS_s_z + SOI_a_r_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(SST_full, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(SST_full))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(SST_full))
qqline(residuals(SST_full), col = "red")Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(SST_full) TDDS_s_z SOI_a_r_z
1.529543 1.529543
# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
TDDS_s_z <- as.numeric(data$TDDS_s_z)
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_SST_full <- Boot(SST_full, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_SST_full) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 1.0852e-16 -0.0016315 0.059488 -0.00096288
TDDS_s_z 6.4909e-02 -0.0003628 0.071978 0.06578146
SOI_a_r_z 1.1184e-01 0.0030682 0.072723 0.11451771
confint(boot_SST_full) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.11907106 0.1159598
TDDS_s_z -0.07670702 0.1985219
SOI_a_r_z -0.03639622 0.2529179
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, TDDS_s_z, SOI_a_r_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
SST_mediate_results <- mediate(model.m = SST_M1,
model.y = SST_full,
treat = "TDDS_s_z",
mediator = "SOI_a_r_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results for the mediation analysis
summary(SST_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.0658 -0.0168 0.14 0.13
ADE 0.0649 -0.0626 0.20 0.39
Total Effect 0.1307 0.0218 0.25 0.03 *
Prop. Mediated 0.5034 -0.3460 2.22 0.16
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs containing zero, neigher sexual disgust nor restricted sociosexual attitudes are significantly positively related to traditionalism.
The coefficients for sexual disgust (.065) and sociosexual attitudes (.112) are quite small.
Mediation Analysis Results:
Even though the proportion mediated is high (50%), it seems that the lack of a relationship between the predictors (over and above one another) and traditionalism means that this large estimated proportion mediated is not reliable.
To attempt to remove the effect of national levels of religiosity and disgust sensitivity from influencing the relationship between disgust and religiosity at the individual level, we will see if we can isolate the individual effect by running a multi-level model with random intercepts. This will account for mean levels of religiosity and isolate individual-level effects from their confounding influence. However, we do not have enough participants in most countries to do this at the country level. Therefore, I will split the nationalities into world region groups. This is admittedly a compromise move, because using national-level groupings would be better. In addition, for many nations, there is not an obvious world region given a lack of representation of participants in countries of those regions. For example, I will exclude participants in areas around the Asia Pacific and participants in South America. We have such nationalities in this sample, but there are too few participants in these regions to justify including them. The grouping that I could come up with that maximizes the number of people in each group while retaining regional coherency was a three-group division: (1) Europe (n = 104), (2) North America (n = 86), and (3) Africa and the Middle East (n = 59). We will split participants up into these regions based on their nationality and fit a multi-level model with a random intercept of region. This will allow us to isolate the individual influence of disgust and religiosity without these three regions’ influence on religiosity scores.
# Load necessary packages
library(lme4) # For fitting the multi-level model
library(sjPlot) # For visualizationLearn more about sjPlot with 'browseVignettes("sjPlot")'.
# Assign regions based on nationality
MLM_data_region <- data %>%
mutate(region = case_when(
nationality %in% c("Canada", "United States", "Mexico") ~ "North America",
nationality %in% c("Belgium", "Austria", "Czech Republic", "Finland", "France",
"Germany", "Greece", "Hungary", "Ireland", "Italy", "Poland",
"Portugal", "Romania", "Slovakia", "Spain", "Sweden",
"United Kingdom") ~ "Europe",
nationality %in% c("Algeria", "Egypt", "Eritrea", "Ghana", "Kenya",
"South Africa", "Saudi Arabia", "Tunisia", "Zimbabwe") ~ "Africa & Middle East",
TRUE ~ NA_character_
))
MLM_data_region <- MLM_data_region %>%
filter(!is.na(region))
# Fit the multilevel model with standardized variables and region as random intercept
MLM_disg_rel_region <- lmer(
CRS_z ~ TDDS_p_z + (1 | region),
data = MLM_data_region,
REML = FALSE
)
# Fixed effects plot
sjPlot::plot_model(MLM_disg_rel_region, type = "pred", terms = "TDDS_p_z")# Random effects plot
sjPlot::plot_model(MLM_disg_rel_region, type = "re")# Faceted scatter plot
ggplot(MLM_data_region, aes(x = TDDS_p_z, y = CRS_z)) +
geom_point(alpha = 0.3, color = "darkblue") +
geom_smooth(method = "lm", se = FALSE, color = "red") +
facet_wrap(~ region) +
theme_minimal() +
labs(title = "Standardized Disgust Sensitivity vs. Religiosity Across Regions",
x = "Standardized Disgust Sensitivity",
y = "Standardized Religiosity")`geom_smooth()` using formula = 'y ~ x'
Looking at the fixed effects plot, we can see that the general tendency is for disgust sensitivity to positively predict religiosity.
Looking at the random intercepts, we can see that while North America and Europe are similar in their levels of religiosity at their mean of disgust sensitivity, Africa and the Middle East is quite high. It is about .75 SD above the mean whereas Europe and North America are about .7 and .3 SD below the mean.
Looking at the faceted scatter plots, it looks like there is a far more distinct positive relationship in Europe—and to a lesser degree in North America—than in Africa and the Middle East between disgust sensitivity and religiosity.
Before moving on to summarizing and testing our model parameter, we will first test some assumptions.
Here we will test the assumptions of homogeneity of variance and normality of residuals for the model by plotting the residuals vs. the fitted values and a histogram and QQ-plot.
# Plotting the residual vs. fitted values
plot(MLM_disg_rel_region, which = 1)# Histogram for normality of residuals
hist(residuals(MLM_disg_rel_region))# QQ-plot for normality of residuals
qqnorm(residuals(MLM_disg_rel_region))
qqline(residuals(MLM_disg_rel_region), col = "red")The residuals vs. fitted values shows two clusters of residuals at the high and low and of the scale, and they are consistent with large differences in religiosity in our three regions. That is, the cluster on the high end of the scale should be Africa and the Middle East. The cluster on the low end of the scale should represent the United States and Europe. This may be a problem for inference, so we will conduct a bootstrapped test of the fixed effect of disgust sensitivity on religiosity.
The histogram and QQ-plot look pretty good.
Now we will conduct a bootstrapped test of the fixed effect of disgust sensitivity on religiosity.
# Set seed for reproducibility
set.seed(1042557)
# Function to extract fixed effects
fixef_fun <- function(model) {
fixef(model)
}
# Perform bootstrapping
boot_results_MLM_disg_rel_region <- bootMer(
MLM_disg_rel_region,
FUN = fixef_fun,
nsim = 1000,
type = "parametric",
cl = cl,
seed = 1042557
)
# Calculate 95% confidence intervals
boot_conf_interval_MLM_disg_rel_region <- apply(boot_results_MLM_disg_rel_region$t, 2, quantile, probs = c(0.025, 0.975), na.rm = TRUE)
# Convert to a more readable format
boot_conf_interval_MLM_disg_rel_region <- as.data.frame(boot_conf_interval_MLM_disg_rel_region)
rownames(boot_conf_interval_MLM_disg_rel_region) <- c("2.5%", "97.5%")
# View the estimates and BCI for fixed effect
summary(MLM_disg_rel_region)Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: CRS_z ~ TDDS_p_z + (1 | region)
Data: MLM_data_region
AIC BIC logLik deviance df.resid
608.1 622.1 -300.0 600.1 245
Scaled residuals:
Min 1Q Median 3Q Max
-2.6269 -0.6843 -0.1220 0.6094 2.7162
Random effects:
Groups Name Variance Std.Dev.
region (Intercept) 0.4390 0.6626
Residual 0.6207 0.7878
Number of obs: 249, groups: region, 3
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.21242 0.38598 0.550
TDDS_p_z 0.09887 0.05160 1.916
Correlation of Fixed Effects:
(Intr)
TDDS_p_z -0.011
print(boot_conf_interval_MLM_disg_rel_region) (Intercept) TDDS_p_z
2.5% -0.5445795 -0.000325435
97.5% 0.9450968 0.199948070
We can see that the random intercept of region accounts for a large amount of variance in religiosity. That is, 44% of the variance in religiosity is attributable to the mean level of religiosity within the region.
However, we can also see that individual disgust sensitivity still is a very marginally significant predictor of religiosity, albeit a smaller size (B = .100, BCI = [.200, -.00003].
Ultimately this analysis suggests that while the relationship between disgust and religiosity may be non-zero, it is being inflated by the presence of Africa and the Middle East in the sample. Similarly, however, the relationship between disgust sensitivy and religiosity appears to be closer to B = 0 within Africa and the Middle East due to restricted range. This suggests that the relationship in Europe and North America may be higher even if accounting for regional variation. To explore the data further, I will fit a similar model, but this time only with North America and Europe as regions.
# Assign regions based on nationality
MLM_data_west <- data %>%
mutate(region = case_when(
nationality %in% c("Canada", "United States", "Mexico") ~ "North America",
nationality %in% c("Belgium", "Austria", "Czech Republic", "Finland", "France",
"Germany", "Greece", "Hungary", "Ireland", "Italy", "Poland",
"Portugal", "Romania", "Slovakia", "Spain", "Sweden",
"United Kingdom") ~ "Europe",
TRUE ~ NA_character_
))
MLM_data_west <- MLM_data_west %>%
filter(!is.na(region))
# Fit the multilevel model with standardized variables and region as random intercept
MLM_disg_rel_west <- lmer(
CRS_z ~ TDDS_p_z + (1 | region),
data = MLM_data_west,
REML = FALSE
)
# Fixed effects plot
sjPlot::plot_model(MLM_disg_rel_west, type = "pred", terms = "TDDS_p_z")# Random effects plot
sjPlot::plot_model(MLM_disg_rel_west, type = "re")# Faceted scatter plot
ggplot(MLM_data_west, aes(x = TDDS_p_z, y = CRS_z)) +
geom_point(alpha = 0.3, color = "darkblue") +
geom_smooth(method = "lm", se = FALSE, color = "red") +
facet_wrap(~ region) +
theme_minimal() +
labs(title = "Standardized Disgust Sensitivity vs. Religiosity Across Regions",
x = "Standardized Disgust Sensitivity",
y = "Standardized Religiosity")`geom_smooth()` using formula = 'y ~ x'
Looking at the fixed effects plot, we can see that the general tendency is for disgust sensitivity to positively predict religiosity.
Looking at the random intercepts, we can see that while North America and Europe are similar in their levels of religiosity at their mean of disgust sensitivity, with North America being slightly more religious.
Looking at the faceted scatter plots,the relationships look similar to the previous plots. There appears to be a positive relationship in both Europe and North America.
Before moving on to summarizing and testing our model parameter, we will first test some assumptions.
Here we will test the assumptions of homogeneity of variance and normality of residuals for the model by plotting the residuals vs. the fitted values and a histogram and QQ-plot.
# Plotting the residual vs. fitted values
plot(MLM_disg_rel_west, which = 1)# Histogram for normality of residuals
hist(residuals(MLM_disg_rel_west))# QQ-plot for normality of residuals
qqnorm(residuals(MLM_disg_rel_west))
qqline(residuals(MLM_disg_rel_west), col = "red")The residuals vs. fitted values doesn’t show much of an issue. There may be more extreme positive residual values and more negative residual values in general. There may also be a slight tendency for residuals to be more negatively extreme towards the positive end of the fitted values. However, there are no glaring patterns indicating a failure of assumptions.
The histogram and QQ-plot look pretty skewed, so we will do a bootstrap inference for the parameter.
Now we will conduct a bootstrapped test of the fixed effect of disgust sensitivity on religiosity.
# Set seed for reproducibility
set.seed(1042557)
# Function to extract fixed effects
fixef_fun <- function(model) {
fixef(model)
}
# Perform bootstrapping
boot_results_MLM_disg_rel_west <- bootMer(
MLM_disg_rel_west,
FUN = fixef_fun,
nsim = 1000,
type = "parametric",
cl = cl,
seed = 1042557
)
# Calculate 95% confidence intervals
boot_conf_interval_MLM_disg_rel_west <- apply(boot_results_MLM_disg_rel_west$t, 2, quantile, probs = c(0.025, 0.975), na.rm = TRUE)
# Convert to a more readable format
boot_conf_interval_MLM_disg_rel_west <- as.data.frame(boot_conf_interval_MLM_disg_rel_west)
rownames(boot_conf_interval_MLM_disg_rel_west) <- c("2.5%", "97.5%")
# View the estimates and BCI for fixed effect
summary(MLM_disg_rel_west)Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: CRS_z ~ TDDS_p_z + (1 | region)
Data: MLM_data_west
AIC BIC logLik deviance df.resid
468.4 481.4 -230.2 460.4 186
Scaled residuals:
Min 1Q Median 3Q Max
-1.4893 -0.6903 -0.2154 0.4676 2.5404
Random effects:
Groups Name Variance Std.Dev.
region (Intercept) 0.008748 0.09353
Residual 0.654947 0.80929
Number of obs: 190, groups: region, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.25509 0.08890 -2.869
TDDS_p_z 0.12206 0.06018 2.028
Correlation of Fixed Effects:
(Intr)
TDDS_p_z 0.090
print(boot_conf_interval_MLM_disg_rel_west) (Intercept) TDDS_p_z
2.5% -0.44074157 0.008464031
97.5% -0.07839034 0.240534554
We can see that the random intercept of region accounts for very little variance in religiosity in this model. Only 0.8% of the variance in religiosity is attributable to the mean level of religiosity within the region.
We can also see that individual disgust sensitivity now a slightly larger predictor of religiosity and is significant (B = .122, BCI = [.241, .008].
This analysis suggests that the relationship between disgust and religiosity in Europe and North America is more robust. This may suggest a relationship between disgust and religiosity that is specific to Christian cultures, because Europe and North America tend to be Christian.
Now we will save the final version of the data frame (data) as a data/final_data.csv.
# Saving the data
write.csv(data, "./data/final_data.csv")